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EXPERIMENTAL  DETERMINATION  OF  DISSIPATED  ENERGY  DENSITY  AS  A 
MEASURE  OF  STRAIN-INDUCED  DAMAGE  IN  COMPOSITES 


Executive  Summary 

The  use  of  composite  materials  in  structural  components  has  increased  dramatically  in  recent  years 
as  their  cost  of  production  continues  to  decline  and  advances  in  composite  design  methodology  become  in¬ 
creasingly  wide  spread.  As  applications  become  more  demanding,  the  need  for  reliable  prediction  of  their 
mechanical  properties  and  behavior  is  becoming  ever  more  important. 

The  basic  obstacle  to  a  comprehensive  understanding  of  failure  behavior  in  composites  is  the  com¬ 
plexity  of  their  observed  mechanical  behavior.  Composites  are  generally  anisotropic,  markedly  nonlinear, 
and,  unlike  metals,  usually  fail  in  an  extremely  complicated,  spatially  diffuse,  noncatastrophic  manner.  The 
nonlinear  mechanical  behavior  of  composites  is  known  to  be  associated  with  damage  accumulation  that 
causes  a  local  decrease  in  stiffness  of  the  material  in  areas  where  damage  is  pronounced.  Moreover,  the  com¬ 
plicated  manufacturing  process  results  in  many  imperfections  that  take  on  a  variety  of  forms  including  fiber 
imperfections,  fiber  misalignment,  and  geometrical  irregularities  in  the  distribution  of  fibers,  voids  or  mi¬ 
crocracks  in  the  matrix  material,  and  the  presence  of  debonded  areas.  These  defects,  combined  with  high- 
stress  fields  near  material  or  geometric  discontinuities,  induce  matrix  cracking,  fiber  breakage,  fiber-matrix 
debonding,  and  delamination,  all  of  which  influence  the  overall  mechanical  properties  of  these  materials. 

Extensive  efforts  to  identify  the  various  modes  of  damage  in  composite  materials  have  been  undertak¬ 
en  in  recent  years.  The  primary  finding  of  most  of  these  investigations  was  that  macroscopic  fracture  was 
usually  preceded  by  an  accumulation  of  the  different  types  of  microscopic  damage  and  occurred  by  the  co¬ 
alescence  of  this  small-scale  damage  into  macroscopic  cracks.  Additionally,  it  was  generally  found  that 
analyses  based  on  classical  fracture  mechanics  did  not  adequately  model  the  damage  effects  and  did  not  pro¬ 
vide  a  satisfactory  degree  of  predictive  capability. 

A  more  practical  approach  to  modeling  failure  behavior  in  composites  is  to  quantify  the  damage  de¬ 
velopment  on  a  continuum  basis  and  relate  it  to  the  material  constitutive  behavior.  The  goal  of  such  an  ap¬ 
proach  is  to  permit  accurate  modeling  of  the  progressive  loss  of  stiffness  and  concomitant  inelastic  behavior 
caused  by  the  underlying  micromechanisms  of  damage.  Such  a  model  for  composite  failure  behavior  re¬ 
quires  a  large  parameter  space  for  its  description  if  it  is  to  represent  physical  fact  in  a  high-fidelity,  objective 
manner.  The  model  must  be  generic  and  phenomenological  in  nature  and  must  rely  on  extensive  sets  of  ex¬ 
perimental  data  to  identify  the  model  parameters.  Only  after  these  parameters  are  determined  can  the  model 
be  related  to  particular  failure  mechanisms. 

The  purpose  of  this  report  is  to  describe  the  initial  stages  of  an  approach  to  modeling  the  mechanical 
behavior  of  composite  materials  consistent  with  the  aforementioned  standards.  This  approach,  which  in  the 
spirit  of  the  prior  discussion  is  largely  based  on  systematic  experimentation,  provides  a  methodology  for 
quantifying  the  failure  behavior  of  a  broad  class  of  composite  materials  and  for  characterizing  the  degree  of 
internal  damage  the  material  undergoes  when  it  is  part  of  a  structural  component  subjected  to  general  load¬ 
ings.  This  approach  also  provides  a  factual  base  against  which  any  other  continuum  damage  approach  can 
be  measured. 

The  procedure  involves  the  determination  of  an  energy  density  dissipation  function  which  has  the  di¬ 
mensions  of  energy  per  unit  volume  and  is  postulated  to  be  a  property  of  the  material.  Its  volume  integral 
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equals  the  energy  dissipated  during  loading  because  of  the  various  internal  failure  events,  and  its  value  at 
any  point  in  the  material  is  regarded  as  a  measure  of  load-induced  internal  damage.  The  energy  dissipation 
function  thus  captures  the  collective  behavior  of  these  failure  mechanisms  without  requiring  an  explicit 
knowledge  of  these  mechanisms,  and,  moreover,  can  also  be  related  to  local  stiffness  changes  that,  as  will 
be  demonstrated  subsequently,  lead  to  a  form  of  nonlinear  structural  behavior.  Section  S  presents  a  detailed 
description  of  the  material  model.  In  our  approach,  the  energy  density  dissipation  function  is  determined 
from  data  obtained  by  means  of  an  extensive  series  of  tests  performed  with  NRL’s  In-Plane  Loader,  which 
is  a  computer-controlled  testing  machine  capable  of  producing  multiple  combinations  of  opening/closing, 
sliding,  and  rotating  boundary  displacements.  A  brief  explanation  of  the  methodology  for  determining 
(identifying)  the  dissipation  function  follows. 

A  representation  of  the  energy  dissipation  function  as  a  function  of  strain  is  chosen  in  terms  of  a  set 
of  known  basis  functions  and  a  set  of  undetermined  coefficients.  The  task  is  then  to  estimate  these  coeffi¬ 
cients.  To  accomplish  this,  a  number  of  test  specimens  are  first  subjected  to  a  range  of  prescribed  combina¬ 
tions  of  in-plane  loads  using  the  In-Plane  Loader.  These  loads  induce  a  wide  range  of  internal  strain  states 
in  the  specimens.  These  states  approximately  encompass  those  expected  in  typical  structural  components 
under  service  loading  conditions.  The  boundary  displacements  and  boundary  loads  are  measured  for  each 
load  combination,  and  are  used  to  compute  the  energy  dissipated  internally  in  the  specimens  caused  by 
strain-induced  damage  for  each  test.  The  strain  fields  associated  with  each  of  the  load  combinations  are  also 
determined  by  means  of  a  linear  finite  element  analysis  using  the  constitutive  properties  of  the  undamaged 
(virgin)  material.  The  computation,  therefore,  neglects  the  effects  of  stiffness  changes  induced  by  the  inter¬ 
nal  damage,  but  should  be  sufficiently  accurate  provided  that  no  large-scale  stress  redistribution  occurs  dur¬ 
ing  the  tests.  Analytic  estimates  of  the  dissipated  energy  are  then  computed  in  terms  of  the  unknown 
coefficients  for  each  test  using  the  representation  of  the  dissipation  function,  the  computed  strain  fields ,  and 
the  fact  that  the  volume  integral  of  the  dissipation  function  equals  the  dissipated  energy.  Finally,  the  unde¬ 
termined  coefficients  in  the  dissipation  function  representation  are  determined  by  minimizing  the  difference 
between  the  experimentally  measured  and  the  analytic  estimates  of  the  dissipated  energy  for  all  the  tests. 
Sections  2  and  3  detail  this  procedure. 

The  aforementioned  procedure  is  a  deconvolution  procedure  in  the  sense  that  the  goal  is  to  extract  a 
material  property  (the  dissipated  energy  density  function)  from  information  obtained  from  the  material  in  a 
particular  physical  configuration,  namely  the  specimen,  and,  thus,  it  is  necessary  to  factor  out  specimen  ge¬ 
ometry  effects.  By  taking  this  approach,  the  resulting  dissipated  energy  density  function  can  be  applied  to 
any  structural  configuration,  and  used  to  develop  material  softening  (stiffness  reduction)  maps  for  any  struc¬ 
ture  (again,  assuming  a  negligible  amount  of  stress  redistribution  caused  by  internal  damage).  Indeed,  in 
Section  4,  we  demonstrate  the  applicability  of  the  model  by  describing  a  computational  approach  packaged 
as  a  material  response  simulator  and  its  use  in  simulating  the  structural  response  of  several  structures  of  in¬ 
terest.  These  are  a  ship’s  mast,  the  test  specimen  itself,  and  a  cylindrical  shell  representing  an  idealized  sec¬ 
tion  of  a  submarine  hull.  A  number  of  different  composite  material  systems  are  considered  and  spatial  maps 
of  the  dissipation  function  (softening  maps)  are  produced  for  various  loading  magnitudes.  These  softening 
maps  illustrate  how  the  energy  consumed  by  the  various  internal  failure  events  is  dissipated  within  the  struc¬ 
ture. 


While  the  energy  dissipation  function  can  be  approximately  determined  without  explicit  knowledge 
of  the  damage-induced  constitutive  nonlinearities  as  discussed  above,  a  more  refined  approximation  re¬ 
quires  a  representation  of  the  full-scale  nonlinear  behavior.  This  information  is  also  of  much  interest  in  its 
own  right  since  it  is  needed  for  stress  analysis  when  analyzing  the  behavior  of  structural  components  loaded 
well  into  the  nonlinear  response  range  where  substantial  load  redistribution  is  expected.  These  issues  are 
addressed  in  Section  5,  where  a  representation  of  the  constitutive  behavior  in  terms  of  the  energy  dissipation 
function  is  developed  along  with  a  refined  scheme  for  identifying  the  dissipation  function  itself. 
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In  summary,  an  approach  is  formulated  to  characterize  failure  behavior  and  degree  of  load  induced 
internal  damage  in  composite  materials  and  structures .  This  approach  is  based  on  a  systematic  experimental 
procedure  to  observe  response  of  composite  materials  subjected  to  multiaxial  load  environment.  The  energy 
dissipated  by  internal  failure  mechanisms  is  used  as  a  measure  of  internal  damage  and  is  characterized  by 
an  energy  dissipation  Junction,  which  is  determined  by  means  of  a  deconvolution  procedure  using  data  pro¬ 
vided  by  NRL’s  automated  In-Plane  Loader.  Use  of  this  information  as  a  failure  analysis  and  prediction  tool 
is  demonstrated  by  simulating  the  structural  response  of  some  naval  structural  components  made  from  sev¬ 
eral  different  composite  materials.  In  addition,  a  general  theory  for  the  derivation  of  the  constitutive  behav¬ 
ior  of  the  damaged  composites  is  presented. 
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EXPERIMENTAL  DETERMINATION  OF  DISSIPATED  ENERGY  DENSITY  AS  A 
MEASURE  OF  STRAIN-INDUCED  DAMAGE  IN  COMPOSITES 


I.  INTRODUCTION 

The  use  of  composite  materials  in  structural  components  has  increased  dramatically  in  recent  years 
as  their  cost  of  production  continues  to  decline  and  advances  in  composite  design  methodology  become  in¬ 
creasingly  wide  spread.  As  applications  become  more  demanding,  the  need  for  reliable  prediction  of  their 
mechanical  properties  and  behavior  is  becoming  ever  more  important. 

The  basic  obstacle  to  a  comprehensive  understanding  of  failure  behavior  in  composites  is  the  sheer 
complexity  of  their  observed  mechanical  behavior.  Composites  are  generally  anisotropic,  markedly  nonlin¬ 
ear,  and,  unlike  metals,  usually  fail  in  an  extremely  complicated  spatially  diffuse  noncatastrophic  manner 
(Badaliance  and  Hill  (1982);  Morris  and  Hetter  (1982);  Masters  and  Reifsnider  (1982)).The  nonlinear  me¬ 
chanical  behavior  of  composites  is  known  to  be  associated  with  damage  accumulation  which  causes  a  local 
decrease  in  stiffness  of  the  material  in  areas  where  damage  is  pronounced.  In  many  fibrous  composite  sys¬ 
tems  damage,  develops  at  applied  loads  much  lower  than  the  design  limit  (Highsmith  and  Reifsnider 
(1982)).  In  materials  that  have  a  laminated  structure,  the  complicated  manufacturing  process  results  in  many 
imperfections  that  take  on  a  variety  of  forms  including  fiber  imperfections,  fiber  misalignment  and  geomet¬ 
rical  irregularities  in  the  distribution  of  fibers,  voids  or  microcracks  in  the  matrix  material,  and  the  presence 
of  debonded  areas.  These  defects,  combined  with  high-stress  fields  near  material  or  geometric  discontinui¬ 
ties,  induce  matrix  cracking,  fiber  breakage,  fiber-matrix  debonding,  and  delamination,  all  of  which  influ¬ 
ence  the  overall  mechanical  properties  of  these  materials.  Extensive  effort  to  identify  the  various  modes  of 
damage  in  composite  materials  has  been  undertaken  in  recent  years  (i.e,.  Aveston  and  Kelly  (1973);  Tsai 
and  Hahn  (1974);  Reifsnider  and  Stinchcomb  (1978);  Reifsnider  (1980);  Highsmith  and  Reifsnider  (1982); 
Dvorak  and  Johnson  (1980)  among  others).  The  primary  outcome  of  these  investigations  was  that  macro¬ 
scopic  fracture  was  usually  preceded  by  an  accumulation  of  the  different  types  of  microscopic  damage  and 
occurred  by  the  coalescence  of  this  small-scale  damage  into  macroscopic  cracks.  Moreover,  it  was  generally 
found  impossible  to  model  the  effects  of  the  damage  using  analyses  based  on  classical  fracture  mechanics 
approaches. 

A  more  practical  approach  to  modeling  failure  behavior  in  composites  is  to  quantify  the  damage  de¬ 
velopment  on  a  continuum  basis  and  relate  it  to  the  material  constitutive  behavior.  The  goal  of  such  an  ap¬ 
proach  is  to  permit  accurate  modeling  of  the  progressive  loss  of  stiffness  and  concomitant  inelastic  behavior 
due  to  the  underlying  micromechanisms  of  damage.  Such  modeling  would  represent  an  extremely  important 
achievement  from  an  engineering  standpoint.  For  example,  consider  a  laminated  composite  plate  with  a  cen¬ 
ter  hole  made  from  a  commercially  available  composite  (i.e.,  AS  1/3501-6  graphite/epoxy)  that  is  loaded  in 
quasi-static  tension.  Damage  patterns  will  develop  around  the  hole  and  near  the  comers  of  the  plate  (Morris 
and  Hetter  (1982))  inducing  degradation  in  material  properties  and  forming  a  precursor  to  general  failure. 
An  examination  of  the  micrographs  readily  indicates  that  a  detailed  micromechanical  analysis  is  out  of  the 
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question  because  the  characteristic  scale  length  of  damage  is  of  the  order  of  |i.m  leading  to  the  conclusion 
that  a  continuum-based  damage  model  would  be  far  more  appropriate  under  these  circumstances. 

The  analysis  of  highly  localized  singular  fields  by  means  of  continuum  approaches  has  a  long  history 
and  there  exist  a  number  of  well-established  techniques  developed  specifically  for  this  purpose.  Most  of 
these  approaches  involve  smearing  out  singular  fields  in  some  fashion  and  replacing  them  by  locally  homo¬ 
geneous  fields,  an  example  being  continuous  dislocation  theory  (i.e.  Bilby,  Bullough  and  Smith  (1955);  Kra- 
tochvil  and  Dillon  (1969),  among  others).  More  generally,  there  have  recently  been  a  plethora  of  continuum 
damage  approaches  based  on  the  idea  of  an  internal  damage  variable  first  introduced  by  L.M.  Kachanov 
(1958)  (e.g.,  Allen,  Harris  and  Groves  (1987),  Bazant  (1986),  L.M.  Kachanov  (1986),  Krajcinovic 
(1984,1986),  Lemaitre  (1987),  Onat  and  Leckie  (1988),  Talreja  (1987)).  A  common  shortcoming  of  most  of 
these  theories  is  that  they  are  usually  not  based  on  observed  facts,  but  rather  tend  to  use  models  defined  over 
low  dimensional  parameter  spaces  that  reflect  preconceived  notions  of  the  dominant  failure  mechanisms. 
Virtually  no  previous  studies  address  this  problem;  indeed,  this  is  one  of  the  primary  purposes  of  this  report. 

An  appropriate  model  for  composite  failure  behavior  requires  a  large  parameter  space  for  its  descrip¬ 
tion  if  it  is  to  represent  physical  fact  in  a  high-fidelity,  objective  manner.  Because  the  phenomena  being  con¬ 
sidered  are  generally  so  complex,  it  seems  futile  to  formulate  models  in  terms  of  the  individual  failure 
events;  rather,  the  chosen  model  must  be  generic  and  phenomenological  in  nature  and  must  rely  on  extensive 
sets  of  experimental  data  to  identify  the  model  parameters.  Only  after  these  parameters  are  determined  can 
the  model  be  related  to  particular  failure  mechanisms. 

This  report  describes  the  initial  stages  of  an  approach  to  modeling  the  mechanical  behavior  of  com¬ 
posite  materials  consistent  with  the  aforementioned  standards.  This  approach,  which,  in  the  spirit  of  the  pri¬ 
or  discussion  is  laigely  based  on  systematic  experimentation,  provides  a  methodology  for  quantifying  the 
failure  behavior  of  a  broad  class  of  composite  materials  and  for  characterizing  the  degree  of  internal  damage 
the  material  undergoes  when  it  is  part  of  a  structural  component  subjected  to  general  loadings.  This  ap¬ 
proach  also  provides  a  factual  base  against  which  any  other  continuum  damage  approach  can  be  measured. 

The  procedure  involves  the  determination  of  an  energy  density  dissipation  Junction  that  has  the  dimen¬ 
sions  of  energy  per  unit  volume  and  is  postulated  to  be  a  property  of  the  material.  Its  volume  integral  equals 
the  energy  dissipated  during  loading  because  of  various  internal  failure  events,  and  its  value  at  any  point  in 
the  material  is  regarded  as  a  measure  of  load-induced  internal  damage.  The  energy  dissipation  function  thus 
captures  the  collective  behavior  of  these  failure  mechanisms  without  requiring  an  explicit  knowledge  of 
these  mechanisms,  and,  moreover,  can  also  be  related  to  local  stiffness  changes  which,  as  will  be  demon¬ 
strated  subsequently,  lead  to  a  form  of  nonlinear  structural  behavior.  In  the  sense  that  this  material  model 
relates  load-induced  internal  damage  to  nonlinear  structural  behavior,  it  is  similar  to  the  models  recently 
used  by  Schapery  (1987, 89a,  89b,  90).  The  present  approach  is  not  limited  to  composites,  but  has  applica¬ 
tion  in  the  general  area  of  the  mechanics  of  irreversible  processes  where  there  is  an  acute  need  for  a  reliable 
and  accurate  procedure  to  experimentally  determine  equivalent  dissipated  energy  density  functions.  This  is¬ 
sue  is  more  fully  addressed  in  several  works  that  deal  with  the  thermomechanical  modelling  of  irreversible 
processes  associated  with  material  behavior  in  general  under  mechanical  loading  (Michopoulos  et  al. 
(1984);  Sih  (1987);  and  Sih  et  al.  (1987)).  Section  5  presents  a  detailed  description  of  the  material  model 
and  its  implementation. 

The  energy  density  dissipation  function  is  determined  from  data  obtained  by  means  of  an  extensive 
series  of  tests  performed  with  NRL’s  In-Plane  Loader  (IPL)  (Mast  et  al.  (1983, 84)).  This  is  a  completely 
computer-controlled  testing  machine  capable  of  producing  multiple  combinations  of  opening/closing,  slid¬ 
ing,  and  rotating  boundary  displacements.  Figure  1  schematically  describes  the  methodology  for  determin¬ 
ing  (identifying)  the  dissipation  function. 
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Fig.  1  —  Block  diagram  of  processes  involved  in  determining  the  Dissipated  Energy  function 


The  left  subpath  in  Fig.  1  describes  the  experimental  part  of  the  identification  process  and  basically 
involves  determining  the  energy  dissipated  in  the  specimens  during  the  IPL  tests.  The  initial  step  consists  of 
using  the  material  and  geometry  specifications  to  manufacture  the  test  specimens.  Next,  the  loading  speci¬ 
fications  are  applied  by  using  the  IPL  to  load  the  specimens  in  the  manner  prescribed.  Finally,  the  measured 
boundary  displacements  and  loads  are  used  to  compute  the  total  energy  dissipated  internally  in  the  speci¬ 
mens  due  to  strain-induced  damage.  Section  2  presents  details  of  this  portion  of  the  identification  process. 
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The  remaining  subpaths  in  Fig.  1  represent  the  numerical  snd  analytical  part  of  the  process.  Initially, 
as  shown  in  the  right  subpath,  the  strain  fields  associated  with  each  of  the  loading  conditions  prescribed  in 
the  loading  specification  are  determined  for  the  given  material  and  specimen  geometry.  This  is  achieved  by 
finite  element  analyses.  A  class  of  functions  with  free  coefficients  is  then  defined  such  that  one  member  of 
this  class  can  represent  an  analytic  representation  of  the  dissipated  energy  as  a  function  of  the  strain  fields. 
Finally,  as  indicated  at  the  bottom  of  Fig.  1 ,  the  particular  function  that  is  best  suited  for  representing  the 
dissipation  energy  is  computed  by  determining  the  values  of  the  free  coefficients  that  minimize  the  differ¬ 
ence  between  the  experimentally  measured  and  the  analytic  representation  of  the  dissipated  energy.  Section 
3  details  this  part  of  the  process. 

The  aforementioned  procedure  is  a  deconvolution  procedure  in  the  sense  that  the  goal  is  to  extract  a 
material  property  (i.e.,  the  dissipated  energy  density  function)  from  information  obtained  from  the  material 
in  a  particular  physical  configuration,  (i.e.,  the  specimen),  and,  thus,  it  is  necessessary  to  factor  out  specimen 
geometry  effects.  By  taking  this  approach,  the  resulting  dissipated  energy  density  function  can  be  applied 
to  any  material  configuration  and  used  to  develop  material  softening  maps  for  any  structure.  In  Section  4, 
we  demonstrate  the  applicability  of  the  model  by  describing  a  computational  approach  packaged  in  a  mate¬ 
rial  response  simulator  and  its  use  in  simulating  the  structural  response  of  a  ship’s  mast  A  number  of  dif¬ 
ferent  composite  material  systems  are  considered  and  spatial  maps  of  the  dissipation  function  are  produced 
for  various  loading  magnitudes.  These  maps  illustrate  how  the  energy  consumed  by  the  various  internal  fail¬ 
ure  events  is  dissipated  within  the  structure  as  a  function  of  load  and  material. 

The  specimen  used  for  the  experimental  measurement  of  the  dissipated  energy  can  also  be  regarded 
as  a  structure  where  the  same  process  can  be  applied.  We  perform  such  an  analysis  in  Section  4.  Although 
this  does  not  seem  to  be  directly  associated  with  naval  structures,  we  chose  to  emphasize  the  fact  that  the 
technique  of  simulating  the  response  of  a  structure  does  not  depend  on  a  particular  shape  or  loading,  but 
rather  on  the  material  response  used  for  its  construction.  At  the  same  time,  this  structure  has  additional  sig¬ 
nificance  because  of  the  essential  role  it  plays  in  the  whole  process.  Finally,  Section  4  also  presents  simula¬ 
tions  of  the  structural  response  of  a  cylindrical  shell  representing  an  idealized  section  of  a  submarine  hull. 

The  procedures  in  our  approach  have  been  implemented  in  a  highly  automated  fashion  by  a  sequence 
of  computational  activities  that  use  the  analyses  described  in  the  following  sections.  Figure  2  shows  a  data 
flow  diagram  view  of  these  processes  and  their  relationships.  The  process  of  acquiring  the  IPL  data  and  the 
computation  of  the  measured  total  dissipated  energy  for  the  test  specimen  is  described  in  Section  2,  and  is 
achieved  through  the  processes  “IPL15”  and  “GETABS,”  the  latter  of  which  evaluates  the  absorbed  and  dis¬ 
sipated  energy  in  the  specimen.  Section  3  describes  the  computation  of  the  dissipated  energy  density  func¬ 
tion.  This  is  performed  by 

•  using  “PATRAN”  to  create  the  complete  finite  element  model  of  the  specimen, 

•  feeding  the  finite  element  data  to  the  “ABAQUS”  finite  element  code  via  the  “PATABA”  data  transla¬ 
tor, 

•  performing  a  linear  elastic  analysis  to  determine  the  strain  field  for  each  point  in  the  specimen,  and 

•  piping  the  measured  dissipated  energy  of  the  specimen  and  the  strain  fields  through  the  “ABAPAT” 
translators  into  “ABSPAT.” 

•  “ABSPAT’  then  computes  the  dissipated  energy  density  function  by  deconvolution  and  uses  an  opti¬ 
mization  scheme  to  minimize  the  difference  between  the  measured  and  the  approximated  values  of 
the  dissipated  energy  as  described  in  Section  3. 

Control  of  the  interprocess  and  intercomputer  communications  is  achieved  through  a  user  interface  module. 
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71  le  simulation  of  material  stiffness  loss  over  any  particular  structure  as  described  in  Section  4  is 
achieved  by 

•  using  “PATRAN”  to  specify  the  finite  element  idealization  and  loading  specification  of  the  structure, 

•  piping  the  data  to  “ABAQUS”  via  the  “PATABA”  data  translator  to  compute  the  strain  fields  associ¬ 
ated  with  the  selected  loading  conditions,  and 

•  piping  the  computed  dissipated  energy  density  function  and  the  strain  field  via  the  “ABAPAT”  trans¬ 
lator  into  “ABSORB”  to  obtain  the  dissipated  energy  distributions  over  the  structure  which  can  then 
be  displayed  using  “PATRAN.” 

Again,  control  of  the  interprocess  and  intercomputer  communications  is  achieved  through  another  user  in¬ 
terface  module. 


Fig.  2  —  Computational  implementation  of  the  current  approach  in  terms  of  the  essential  modules  and  their 
functional  relationship 
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More  detailed  versions  of  the  computational  descriptions  pertaining  to  all  the  aforementioned  pro¬ 
cesses  are  presented  in  the  corresponding  sections  (2, 3  and  4)  of  the  report. 

While  the  energy  dissipation  function  can  be  approximately  determined  without  explicit  knowledge 
of  the  damage-induced  constitutive  nonlinearities  as  discussed  above,  a  more  refined  approximation  requires 
a  representation  of  the  full-scale  nonlinear  behavior.  This  information  is  also  of  much  interest  in  its  own 
right  since  it  is  needed  for  stress  analysis  when  analyzing  the  behavior  of  structural  components  loaded  well 
into  the  nonlinear  response  range  where  substantial  load  redistribution  is  expected.  Section  5  addresses  and 
develops  a  representation  of  the  constitutive  behavior  in  terms  of  the  energy  dissipation  function  along  with 
a  refined  scheme  for  identifying  the  dissipation  function  itself. 

2.  DATA  ACQUISITION  THROUGH  THE  IN-PLANE  LOADER  SYSTEM 
2.1  Material  Specification 

Composites  associated  with  naval  applications  range  through  a  wide  variety  of  materials.  Each  differ¬ 
ent  combination  of  matrix,  fiber,  fiber  coating  (for  matrix-fiber  interface),  layup  angle,  stacking  sequence, 
etc.,  corresponds  to  a  different  material.  Our  approach  is  specifically  tailored  to  organic  matrix  composites. 
In  general,  the  particular  material  specification  for  an  application  in  a  structure  depends  not  only  on  the  me¬ 
chanical  considerations  but  also  on  a  host  of  additional  considerat;ons  such  as  cost,  electromagnetic  and 
thermal  properties.  Approximately  100  material  systems  with  graphite  fibers  ranging  from  AS1  to  IM7 
(Hercules  Corp.)  and  several  thermoset  and  thermoplastic  organic  polymers  have  been  tested  and  character¬ 
ized  with  the  present  approach.  Here  we  limit  our  description  to  the  materials  selected  for  the  applications 
discussed  in  this  report.  A  brief  description  of  these  materials  is  given  in  Table  1  in  terms  of  the  matrix,  the 
fiber,  and  the  angle  of  layup.  Appendix  A  is  a  complete  list  of  all  the  materials  tested  up  until  now. 

Most  of  the  materials  were  received  in  the  form  of  panels  out  of  which  the  specimens  were  cut  and 
machined  to  the  specimen  specification. 

Table  1  —  Fragment  of  the  Materials  Database  Describing  Four  of  the  101  Different 
Materials  Used  in  the  Present  Study 


ITEM 

RESIN 

FIBER 

LAYUP 

SUPPLIER 

FILE 

PROJECT 

001 

3501-6 

AS1 

+/-  15 

Hercules 

mt-1 

NAV-SEA 

002 

. 

• 

+/-30 

. 

mt-2 

. 

007 

. 

AS4 

+/-60 

. 

mt-7 

. 

040 

PEEK 

. 

+/-  60 

IC1 

mt-40 

. 

2 J2  Specimen  Specification 

The  specimen  geometry  shown  in  Fig.  3  was  designed  to  satisfy  the  following  requirements: 

•  The  characteristic  dimensions  be  large  enough  relative  to  fiber  diameter  and  lamina  thickness  to  en¬ 
sure  that  the  material  could  be  analyzed  as  either  a  single  mechanically  equivalent  homogeneous 
anisotropic  monolithic  material,  or  a  collection  of  layers  of  varying  orientations  of  such  materials. 

•  The  overall  specimen  size  be  small  enough  to  keep  material  costs  at  a  manageable  level. 
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•  A  strain  riser  be  present  to  guarantee  that  high-strain  regions  occur  well  away  from  ail  specimen 
boundaries. 

Figure  3  shows  the  single  edge-notched  specimen  that  resulted.  It  is  important  to  note  that  the  primary 
function  of  the  notch  is  to  act  as  a  strain  riser  to  ensure  satisfaction  of  the  last  requirement  above,  and  not  to 
mimic  a  crack  in  the  fracture  mechanics  sense. 


Fig.  3  —  (a)  Schematic  of  the  specimen  and  the  associated  areas,  (b)  photos  of  representative  specimen  before  testing 
(left)  and  after  testing  (right)  for  material  2  with  3501-6  resin  AS1  fiber  and  +/-  30°  layup  and  loading  path  1 1 


2.3  Loading  Specification 

The  experimental  system  used  for  measuring  th'  dissipated  energy  in  the  test  specimens  is  the  In- 
Plane  Loader  System  (IPLS).  Appendix  B  presents  a  defied  description  of  this  testing  system.  The  primary 
component  of  the  IPLS  is  a  custom-made,  fully  automated  testing  machine  (the  In-Plane  Loader  or  IPL) 
shown  in  Fig.  4.  The  remaining  components  of  the  IPLS  are  a  computer  that  fully  controls  the  IPL  and  a 
graphics  processor  that  is  used  for  various  postprocessing  operations  in  conjunction  with  the  computer  sys¬ 
tem. 


The  objective  of  the  IPLS  is  to  control  the  rigid  body  motion  of  the  boundary  of  the  specimen  that  is 
held  by  the  movable  grip.  Because  the  actuators  are  constrained  to  move  in  a  plane  parallel  to  the  specimen, 
the  resulting  motion  involves  only  three  degrees  of  freedom  relative  to  any  frame  of  reference  on  that  plane. 

The  relation  between  the  prescribed  actuator  displacements  and  the  resulting  grip  motion  is  illustrated 
in  Fig.  B3(a)  and  B3(b)  of  Appendix  B.  It  is  important  to  note  here  is  that  the  grip  motion  can  be  resolved 
into  three  basic  components:  sliding  ux,  opening/closing  uy,  and  rotation  CD.  Specified  combinations  of  ac¬ 
tuator  displacements,  therefore,  map  into  particular  combinations  of  these  three  basic  motions. 

It  is  convenient  to  use  a  reference  frame  located  at  the  initial  position  of  the  notch  tip  for  both  the 
boundary  displacements  (denoted  by  u  relative  to  this  frame)  and  the  resulting  reaction  forces  £.  Figure 
B3(c)  of  Appendix  B  illustrates  these  displacements  and  tractions  and  the  resulting  deformation  of  the  spec- 
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Fig.  4  —  (a)  Side  view  of  the  In-Plane  Loader,  (b)  View  of  the  grip  area  of  the  In-Plane  Loader 

The  decomposition  of  the  applied  displacements  relative  to  this  frame  into  a  sliding  motion  Uq,  an 
opening/closing  motion  u{,  and  a  rotational  motion  u2  is  shown  in  Fig.  B3  (d,e,f)  of  Appendix  B  and  Fig. 
5(b,  c  and  d),  respectively.  For  dimensional  homogeneity,  u2  is  defined  as  the  length  of  the  arc  traveled  by  a 
point  1  inch  away  from  the  notch  tip  rigidly  connected  with  the  moving  grip  along  the  direction  of  the  rota¬ 
tion,  instead  of  using  the  actual  rotation  in  radians. 
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Fig.  5  —  View  area  of  the  1PL  specimen  in  the  grips  (a)  before  applying  any  loading,  (b)  after  applying  some  translation 
of  the  upper  grip  parallel  to  the  x-  axis,  (c)  parallel  to  the  y-  axis,  and  (d)  after  applying  rotation  about  the  origin  of  the  x- 
y  frame  of  reference 


Subsequent  analysis  requires  that  we  compute  the  energy  dissipated  within  the  material  at  a  discrete 
set  of  observation  points  as  the  specimen  is  loaded  by  applying  a  predetermined  series  of  boundary  displace¬ 
ments  if.  Here  p=\,...ji  and  n  denotes  the  number  of  observation  points.  The  details  of  these  computations 
will  be  discussed  in  the  next  section;  here  we  attend  to  the  sampling  strategy.  It  is  advantageous  to  think  in 
terms  of  a  three-dimensional  displacement  space  with  coordinates  (uq,  uj,  u2).  The  issue  then  is  how  to  se¬ 
lect  a  representative  family  of  paths  that  cover  the  space  and  how  to  sample  along  each  path. 
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It  is  not  expected  to  observe  any  significant  path-dependent  behavior  during  this  initial  loading 
phase.  By  this  we  mean  that  we  expect  the  mechanical  response  of  the  material  at  any  point  in  this  space  to 
depend  only  on  the  current  state  of  internal  strain  and  to  be  independent  of  the  particular  path  followed  to 
achieve  this  state.  This  type  of  path  independence  is  consistent  with  the  observations  of  Schapery  and  Lam- 
bom  (1988)  and  greatly  simplifies  matters  because  it  permits  us  to  cover  the  displacement  space  with  a  fam¬ 
ily  of  loading  paths  selected  solely  on  the  basis  of  convenience. 

Towards  this  end,  it  was  decided  to  cover  the  boundary  displacement  space  with  a  set  of  15  uniformly 
distributed  radial  loading  paths  as  indicated  in  Fig.  6.  Note  that  because  of  geometry  and  material  symmetry 
about  the  x  axis  (Fig.  5),  only  the  half-space  corresponding  to  positive  sliding  displacement  (uq>0)  need  be 
considered.  The  required  set  of  observation  points  is  generated  by  sampling  along  each  path  at  50  distinct 
points  starting  from  0  mils  and  terminating  at  a  maximum  of  50  mils  of  displacement  yielding  a  total  of  750 
points  per  material  system.  Then  a  particular  test  in  that  the  actuator  motions  are  continuously  varied  cor¬ 
responds  to  a  specific  path  in  this  space.  This  path  can  be  represented  by  a  vector  originating  from  the  origin 
of  the  space  and  components  given  by  u  •  =  rai  for  i  =  0, 1,  2  ,  where  a(- are  the  coordinates  of  the 
unit  vectors  along  the  loading  path  direction,  and  r  is  a  scalar  multiple  denoting  the  proportionality  of  the 
path  and  ranging  from  0  to  0.05  inches  in  steps  of  0.001  inches  corresponding  to  the  successive  observation 
points.  Only  15  specimens  are  required,  and  50  observations  per  loading  path  are  obtained  from  a  single 
specimen.  Table  2  shows  the  coordinates  ai  of  the  unit  vectors  for  the  selected  loading  paths. 


Fig.  6  —  Definition  of  the  proportional  loading  paths  in  the  boundary  displacement  space  (loading  path  1 1 
is  a  representative  case),  and  the  uniform  distribution  of  the  IS  paths  used  in  the  present  methodology 
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Table  2  —  Coordinates  of  Load  Path  Unit  Vectors  for  each  one 
of  the  15  Loading  Cases  Applied  by  the  IPL 


LOADING 

<*0 

<*1 

a2 

1 

0.000 

-0.707 

-0.707 

2 

0.000 

-1.000 

0.000 

3 

0.000 

-0.707 

0.707 

4 

0.500 

-0.500 

-0.707 

5 

0.707 

-0.707 

0.000 

6 

0.500 

-0.500 

0.707 

7 

0.707 

0.000 

-0.707 

8 

1.000 

0.000 

0.000 

9 

0.707 

0.000 

0.707 

10 

0.500 

0.500 

-0.707 

11 

0.707 

0.707 

0.000 

12 

0.500 

0.500 

0.707 

13 

0.000 

0.707 

-0.707 

14 

0.000 

1.000 

0.000 

15 

0.000 

0.707 

0.707 

The  locus  of  the  end  points  of  all  loading  paths  for  the  same  increment  is  a  half-sphere  as  shown  in 
Fig.  6,  where  loading  path  1 1  at  an  arbitrary  increment  is  presented  as  an  example. 

2.4  Testing  Procedure 

The  IPL  is  instructed  to  apply  loads  along  the  designated  loading  path  for  each  specimen.  Because 
each  test  is  performed  twice  to  establish  the  degree  of  reproducibility,  30  (2  per  loading  path)  specimens  are 
required  for  analyzing  a  single  material.  After  the  IPL  calibration  phase  is  complete,  a  stack  of  thirty  spec¬ 
imens  from  the  same  material  is  loaded  on  the  magazine  of  the  automated  specimen  feeder  of  the  IPL  and 
the  testing  phase  follows.  The  following  sequence  of  events  is  repeated  for  each  specimen. 

•  The  specimen  is  fed  from  the  specimen  magazine  into  the  grip  area  of  the  IPL. 

•  Pressure  is  applied  on  the  hydraulic  grips  to  clamp  the  specimen  with  the  appropriate  force  to  prevent 
slippage  and  crushing  of  the  specimen. 

•  The  operator  uses  the  digital  imaging  system  of  the  IPL  to  control  the  positioning  of  a  crosshair  on  the 
image  of  the  area  between  the  grips  as  acquired  by  the  overhead  video  camera.  Notch  tip  location  and 
selected  points  between  the  grips  are  thus  automatically  fed  into  the  computer  of  the  IPLS,  to  establish 
the  frame  of  reference. 

•  The  remainder  of  the  test  is  performed  automatically,  with  the  computer  control  system  measuring 
specimen  thickness  and  controlling  the  grip  pressure,  selecting  and  applying  incrementally  the  load¬ 
ing  path.  The  system  monitors  and  stores  the  boundary  forces  and  displacements  for  each  loading  step, 
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and  using  this  data  the  total  energy  dissipated  in  the  specimen  for  each  loading  increment  in  the  load¬ 
ing  path  is  computed  and  stored  along  with  the  force  displacement  data. 

2.5  Computational  and  Data  Description 

Figure  7  is  a  data  flow  diagram  showing  the  computational  process  corresponding  to  the  left  subpath 
of  Fig.  1 .  The  data  from  the  A/D  and  D/A  converters  are  directed  to  module  “IPL15”  through  the  Input/Out- 


Fig.  7  —  Block  diagram  of  the  computational  modules  and  resources  involved  with  the  determination  of  the 
measured  total  dissipated  energy 


put  driver  module.  “IPL15”  controls  the  data  acquisition  and  control  process  by  controlling  the  following 
operations: 

•  it  controls  the  positioning  of  the  actuators  and  the  associated  hardware; 

•  it  collects  the  force  and  displacements  (f, «)  at  the  50  sampling  points  along  the  loading  path  from  the 
transducers; 

•  it  then  transforms  them  into  those  referring  to  the  frame  of  reference  attached  on  the  tip  of  the  notch 
as  in  Fig.  5(a); 

•  it  computes  and  displays  the  total  and  instantaneous  values  of  the  total  dissipated  energy  lost  in  the 
specimen  due  to  damage  according  to  the  methodology  described  in  Section  2.6; 

•  it  finally  stores  all  measured  quantities  in  files  (one  per  material).  Here  the  term  “material”  is  used  for 
a  laminate  of  specific  fiber,  matrix,  and  lamination  angle. 

Each  file  contains  the  data  for  all  specimens  and  loading  paths  for  a  given  material.  Each  file  consists 
of  a  header  and  a  series  of  blocks  that  contain  information  pertaining  to  each  specimen.  The  header  is  at  the 
beginning  of  the  file  and  contains  the  original  filename.  Each  specimen  block  has  its  own  header  (which  is 
seven  records  long)  and  fifty  data  records  (one  for  each  experimental  point  on  the  loading  path).  The  header 
for  each  specimen  contains  all  the  pertinent  information  associated  with  the  corresponding  specimen  and 
the  respective  loading  path.  Figure  8  presents  a  section  of  one  of  these  files  that  corresponds  to  the  second 
material  of  Table  1 ,  i.e.,  graphite  epoxy  with  AS  1  fiber  and  3501-6  resin  with  layup  angle  +/-  30  °.  The  fields 
in  the  data  records,  from  left  to  right,  are  experimental  point  number,  total  dissipated  energy  D,  boundary 
displacements  (uq,  u\,  and  u^),  and  boundary  forces  (ig,  t j,  and  t^. 


12 


MASTETAL. 


Fig.  8  —  Typical  parts  of  file  containing  the  measured  through  the  In-Plane  Loader  data 

Another  process  called  “GETABS”  extracts  the  dissipated  energy  values  from  the  aforementioned  file 
and  stores  them  into  files  ready  to  be  used  from  the  resources  available  for  the  determination  of  the  dissipat¬ 
ed  energy  density  function  as  described  in  Section  4. 
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All  these  processes  were  originally  executed  on  a  PRIME  750  mv  computer,  but  recently  an  effort 
has  been  undertaken  to  recast  them  in  a  more  contemporary  form  on  the  NeXT  cube  computing  platform. 

2.6  Computation  of  the  Dissipated  Energy 

In  this  section  we  discuss  in  detail  the  procedure  for  computing  the  energy  dissipated  internally  in  the 
specimen  during  an  IPL  test.  We  employ  the  following  convention  in  the  rest  of  this  report:  superscript  p  (p 
=  1,...,50)  denotes  an  experimental  point  on  a  given  loading  path;  1(1  =  1,...,15)  denotes  a  particular  loading 
path;  r  =  1,2  denotes  a  particular  test  repetition;  and  subscript  i  (i  =  0,1,2)  denotes  the  components  of  the 
boundary  forces  and  displacements. 

Consider  first  the  total  energy  imparted  to  the  specimen  for  each  experimental  point  p  as  function  of 
the  boundary  displacement  vector  u.  These  energies,  which  we  denote  as  W*\  can  be  computed  from  the 
measured  boundary  forces  and  displacements  by  using  the  standard  definition  of  mechanical  work.  Thus 

/ 

=  W^Qr)  =  jf-ds, 

0 

or  in  component  form. 

Hi  uS  Ut  Hi 

Wp=  /<?*,=  jlSdl0+  /»?*,+  jt£ds2 
0  0  0  0 

Each  term  in  Eq.  (2)  represents  the  energy  stored  in  the  specimen  due  to  one  of  the  measured  force-dis¬ 
placement  pairs  (tf,  uP)  and  is  denoted  by  W?.  Hence,  we  may  express  W**  as  the  sum 

Vf  =  W^+WP+W?.  (3) 


The  integrals  in  Eq.  (2)  can  be  approximated  using  a  discrete  representation,  i.e., 

Wg  =  \  X  <ro  +*o"1)  («o““o_1) 

a  =  1 

W?  =  5  £  (I (”|+<T ')  -  «‘r ' )  (4) 


W?  =  jt  <»I  +*5~')  <“2  -»v'  >. 

a  =  1 


(1) 


(2) 


or  in  general 
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!  P 

wf-  «*■  *  =  0,1,2.  (5) 

a  =  1 

This  approximation  process  is  graphically  depicted  in  Fig.  9  where  the  quantities  in  the  summation 
are  defined  geometrically  on  a  representative  traction  vs  displacement  graph. 

The  area  below  the  curve  represents  the  energy  W?  for  loading  up  to  point  p  as  the  sum  of  all  trape¬ 
zoids  of  area  *  (m*  _w*- 1  )  • 

The  total  energy  W*’ can  now  be  computed  by  substituting  Eq.  (5)  into  Eq.  (3)  and  taking  the  sum  of 
the  three  terms.  Thus, 


for?  =1,2 . 50.  (6) 

«  =  0a  =  1 

Consider  now  the  situation  when  the  specimen  unloads.  We  assume  that  when  the  specimen  is  unload¬ 
ed  proportionally  from  any  point  a  of  the  curve  in  Fig.  9,  the  unloading  occurs  elastically  along  a  line  from 
the  point  a  to  the  origin.  This  assumption  is  largely  based  on  observations  made  during  the  early  years  of 
1PL  use  when  unloading  data  were  collected;  it  was  found  then  that  almost  all  specimens  tested  unloaded 
linearly  and  showed  a  maximum  2%  permanent  displacement  set.  These  permanent  displacements  are  likely 
a  consequence  of  initial  stresses  introduced  during  manufacture  and  were  deemed  small  enough  to  ignore. 
This  characteristic  non-metallike  behavior  during  unloading  has  also  been  observed  by  other  investigators  ( 
Shapery  and  Lam  bom  (1988))  and  it  is  currently  attributed  to  the  propensity  of  organic  composites  to  form 
microcracks  instead  of  flowing  plastically. 


Fig.  9  —  Schematic  representation  of  the  integration  scheme  used  for  the  compute 
W?  as  the  sum  of  the  rectangular  slices  under  the  load-displacement  curve  and  of 
the  energies  associated  with  the  loading  of  the  specimen  up  to  the  point  p 
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During  unloading,  a  portion  of  the  total  energy  W?  is  recovered  because  of  the  elastic  nature  of  the 
unloading  response.  This  recoverable  energy,  denoted  by  a  componentwise  basis  as  fl?,  is  simply  the  trian¬ 
gular  area  under  the  line  from  a  to  the  origin,  or 

for  p  =  1,2,...,  50.  (7) 

Knowing  the  available  energies  W?  and  the  recoverable  energies  flf.it  is  a  simple  matter  to  compute 
the  energy  D p  dissipated  in  the  specimen.  Thus,  on  a  componentwise  basis,  the  energy  dissipated  in  the 
specimen  as  it  is  loaded  up  to  point  p  is  simply  the  difference 


Df=Bf-Af,  for  p=  1,2 . 50,  (8) 

and,  in  graphical  terms,  is  just  the  shaded  area  under  the  load-displacement  line  in  Fig.  9.  Substituting  the 
expressions  for  the  total  and  available  energies  from  Eqs.  (5)  and  (7)  for  each  component  i  in  the  relation 
above,  therefore,  gives  for  Df 


*r- Jif'.+o  -J«.  for  p~ . 5°- 


m 


a  =  1 

The  total  dissipated  energy  D?  (the  energy  which  appears  in  the  second  column  of  the  acquired  data  file 
shown  in  Fig.  8)  can  now  be  computed  by  summing  Z>f  componentwise  and  results  in 


o' -if  i (/■,+if-1) («•-«**') for  P-i,* . 


50. 


(10) 


i  =  0a  =  1 


i  =  0 


The  dissipated  energy  computed  this  way  can  be  considered  to  be  a  measured  value  since  it  is  derived 
directly  from  measured  quantities  and  the  only  sources  of  error  are  from  the  discrete  numerical  integration 
described  and  the  quantization  error  of  the  data  acquisition  process.  However,  the  very  operation  of  integra¬ 
tion  also  serves  as  a  noise/error  reduction  methodology  and  works  more  in  a  beneficial  rather  than  a  detri¬ 
mental  manner  in  this  case.  As  far  as  we  know,  these  energy  values  represent  the  only  attempt  to  measure 
dissipated  energy  caused  by  mechanical  deformation  in  composite  materials. 


Figure  10  is  a  graphical  representation  of  the  boundary  displacements  and  tractions  as  well  as  the  as¬ 
sociated  total  and  dissipated  energies  for  the  first  specimen  for  loading  path  1 1  for  material  AS  1-30  (second 
row  of  Table  1).  Loading  path  1 1  does  not  have  any  rotation  applied  on  the  specimen  as  can  be  seen  from 
Table  2  and  from  Fig.  6.  This  results  in  the  apparent  zeroing  of  the  curve  in  Fig.  10(a.2),  which  consequently 
results  in  the  blank  energy  graphs  10(c.2)  and  10(d.2).  However,  due  to  the  other  two  displacement  compo¬ 
nents,  there  is  a  force  reaction  in  the  rotation  direction  as  is  shown  in  Fig.  10(b.2). 


The  dissipated  energy  curves  start  with  values  very  close  to  zero  and  suddenly  rise  steeply.  This  hap¬ 
pens  as  the  corresponding  displacements  exhibit  an  abrupt  increase  and  the  tractions  reach  a  maximum  value 
and  then  rapidly  decrease.  This  point  can  be  interpreted  as  defining  a  transition  between  purely  linear  elastic 
material  behavior  (i.e.,  no  energy  dissipation)  and  the  point  at  which  the  material  first  exhibits  nonlinear  be¬ 
havior  with  accompanying  energy  dissipation. 

Summation  of  the  energies  for  the  three  components  in  Figs.  10(c.0),  10(c.l),  and  10(c.2)  yields  Fig. 
10(e)  for  the  energy  absorbed  from  all  three  component  combinations  in  the  specimen.  Similarly,  summation 
of  the  dissipated  energy  components  in  Figs.  10(d.0),  10(d.l),  and  10(d.2)  yields  the  total  energy  lost  in  the 
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Fig.  10  —  Evolution  of  displacements  (a.0,  a.1,  a.2),  tractions  (b.0,  b.l,  b.2),  total  absorbed  energies  (c.0,  c.  1,  c.2),  dissipated 
energies  (d.0,  d.l,  d.2),  total  absorbed  energy  (e),  and  total  dissipated  energy  (f),  vs  the  magnitude  of  the  experimentally  im¬ 
posed  displacement  vector  Hull  of  loading  path  1 1  for  the  first  specimen  of  material  2  ( AS1, 3501-6,  +/-  30  0 ) 
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specimen  shown  in  Fig.  10(f).  The  morphology  of  the  displacement,  traction,  and  energy  evolution  curves 
for  the  rest  of  the  loading  paths  is  not  necessarily  identical  to  that  of  loading  path  11. 

2.7  Statistics  of  the  Process 

To  compare  the  repeatability  of  the  process.  Fig.  1 1  shows  the  dissipated  energy  evolutions  for  both 
specimens  used  on  loading  path  1 1  together  with  their  average  and  absolute  difference  evolutions. 

As  it  seen  in  Fig.  11(d)  the  specimens  behave  almost  identically.  A  small  difference  of  the  order  of 
12%  appears  at  the  last  experimental  point,  very  deep  in  the  nonlinear  region  where  the  enormous  amount 
of  damage  has  made  the  material  more  unstable  than  it  is  in  the  linear  region.  The  root  mean  square  error 
(RMS)  between  the  two  values  of  the  dissipated  energy  for  the  two  specimens  from  the  same  material  has 
never  been  found  to  exceed  5%  for  all  materials  tested. 

An  entire  IPL  test  takes  about  10  seconds,  five  of  which  are  spent  installing  the  specimen  in  the  grips. 
As  a  result,  the  specimen  testing  rate  can  be  360  specimens/hour.  At  this  rate,  data  can  be  collected  for  12 
different  materials  per  hour.  The  overall  daily  production  rate  for  an  8-hour  day  is,  therefore,  960  specimens, 
or  96  different  materials,  or  24  materials  systems  (since  we  use  4  layup  angle  combinations  for  each  fiber- 
resin  combination).  The  total  number  of  experimental  points  per  fiber-resin  combination  is  6000.  Each  load¬ 
ing  path  corresponds  to  2.4  KB  of  data  while  288  KB  are  acquired  for  each  material  system.  The  daily 
throughput  capacity  of  the  acquisition  process  approaches  20.74  MB/day. 


Fig.  1 1  —  Dissipated  energies  for  specimen  1  (a),  specimen  2  (b),  their  average  value  (c),  and  the  absolute  difference  between 
the  two  (d)  for  material  2  (  AS  1/6501-6  [+/-  30  °  ])  and  loading  path  1 1 
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3.  DISSIPATED  ENERGY  DENSITY  FUNCTION  DETERMINATION 

This  section  is  primarily  concerned  with  the  right  and  the  bottom  subpaths  of  Fig.  1.  These  subpaths 
represent  the  flow  of  activities  related  to  the  determination  of  the  dissipated  energy  density  function  <j>  .  It 
is  important,  however,  to  clarify  the  range  of  hypotheses  and  assumptions  employed  in  the  present  analysis 
to  illustrate  how  our  procedure  relates  to  the  so-called  inverse  approach  in  mechanics  of  materials  research. 


3.1  Computational  Procedures 

Figure  12  describes  the  computational  process  used  to  obtain  the  dissipated  energy  density  function 
in  terms  of  its  evaluation  on  discrete  points  over  the  strain  space. 


Fig.  12  —  Block  diagram  of  the  modules  used  to  compute  the  dissipated  energy  coefficients  for  each  material,  along  with 
the  computational  resources  utilization 


A  user  interface  that  controls  all  interprocess  and  intercomputer  communication  is  the  mechanism  by 
which  the  user  controls  the  overall  process.  The  “PATRAN”  solid  modeling  package  is  used  to  define  the 
geometry,  material,  and  loading  specifications  for  the  finite  element  model  of  the  specimen  structure,  which 
in  turn  is  fed  into  the  “ABAQUS”  finite  element  code  via  the  “PATABA”  format  translator.  The  output  from 
the  linear  strain  analysis  from  “ABAQUS”  is  then  translated  from  “ABAPAT”  into  files  containing  die  strain 
and  displacement  distributions  for  the  corresponding  loading  cases  that  are  displayable  by  “PATRAN.”  The 
“ABSPAT”  module  developed  at  NRL  then 

•  uses  the  analysis  results  for  all  the  loading  cases  as  well  as  the  experimentally  measured  values  of  the 
dissipated  energy  from  the  files  discussed  in  Section  2.3  to  synthesize  the  quantities  necessary  for  es¬ 
timating  the  125  coefficients  that  describe  the  dissipated  energy  function, 

•  carry  out  the  estimation,  and 

•  generate  one  file  per  material  system  that  contains  the  125  coefficients. 
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3 2  Strain  Field  Determination 

The  measured  absorbed  energy  is  a  function  of  specimen  geometry,  loading  condition,  and  material. 
To  factor  out  the  influence  of  specimen  geometry  and  loading  and  obtain  a  dissipated  energy  function  that 
depends  solely  on  material,  it  is  necessary  to  develop  a  deconvolution  procedure.  A  description  of  our  ap¬ 
proach  to  achieving  this  goal  follows. 

For  every  structure  of  a  certain  geometry  s,  consisting  of  a  particular  material  m,  and  under  a  specific 
loading  condition  /,  a  strain  state  exists  at  any  point  described  by  the  position  vector  x.  The  strain  vector 
corresponding  to  this  state  with  the  classical  in-plane  components  (£  =  (£xx,  £  ,  £  )  )  is  expected  to  be 
a  function  of  s,  m,  /,  and  x;  that  is, 


£  =  e(s,  l,  m,x) .  (11) 

When  the  structure  under  consideration  is  a  specific  one,  such  as  the  specimen  shown  in  Fig.  2,  a  sub¬ 
script  0  for  both  the  corresponding  geometry  and  material  is  used  to  refer  to  this  structure  (sg)  and  material 
(mg)  thereafter.  It  will  become  evident  from  the  analysis  that  the  strain  field  for  each  one  of  the  loading  cases 
needs  to  be  determined.  In  the  ideal  case,  these  fields  should  be  experimentally  measured  along  with  bound¬ 
ary  tractions  and  displacements  in  the  IPL  tests.  This  will  be  a  feature  of  the  new  generation  of  NRL’s  auto¬ 
mated  loading  machines,  however,  the  current  generation  has  no  such  capabilities.  The  results  and  analysis 
presented  here  are  based  on  an  alternative  approach.  The  approach  used  in  this  study  to  estimate  the  strain 
field  was  to  discretize  the  geometry  of  the  specimen  into  150  elements  as  shown  in  Fig.  13,  systematically 
apply  the  loading  and  boundary  conditions  applied  by  the  IPL  on  the  specimen,  and  compute  the  in-plane 
strains  by  running  an  appropriate  number  of  linear  elastic  finite  element  analyses  using  the  code 
"ABAQUS.”  The  selection  of  the  linear  analysis  was  based  on  the  fact  that  in  very  few  of  the  experiments 
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were  specimens  observed  to  fail  in  a  way  that  would  indicate  that  very  large  amounts  of  energy  had  been 
dissipated  by  internal  material  damage.  In  addition,  our  initial  purpose  was  to  establish  an  analytical  tool 
based  on  experimental  evidence  for  predicting  the  onset  of  nonlinear  response  and  an  approximation  of  this 
response,  rather  than  the  exact  nonlinear  material  behavior.  The  three  strain  components  corresponding  to 
the  centroids  of  the  elements  of  the  specimen  structure  sq  were  stored  in  files  organized  according  to  loading 
case.  For  practical  purposes,  instead  of  running  15  loading  cases  (one  for  each  loading  path),  only  the  cases 
corresponding  to  the  loading  paths  ( 1 , 0, 0),  (0, 1 , 0),  and  (0, 0, 1 )  corresponding  to  pure  shear,  opening,  and 
rotation  displacements,  respectively,  were  run.  This  was  because  all  other  cases  could  be  synthesized  by  lin¬ 
ear  superposition  according  to  Table  2. 

Figure  13  shows  the  finite  element  mesh  idealization  of  the  specimen.  Although  initially  the  entire 
area  between  the  grips  had  be  discretized,  subsequent  analysis  indicated  that  the  area  close  to  the  left  side 
of  the  specimen  (the  side  of  the  notch)  remains  unstressed  for  virtually  all  loading  cases.  Thus,  only  the  view 
area  (as  shown  in  Fig.  2)  was  actually  analyzed. 

3.3  Dissipated  Energy  Density  Function 

We  assume  that  there  exists  a  scalar  function  $  expressing  a  measure  of  the  dissipated  energy  density 
per  unit  of  volume  of  the  material,  which  only  depends  on  the  strain  vector  £  and  the  material  used  in  the 
structure,  i.e.,  <J>  (e,  m). 

We  furthermore  assume  that  this  function  can  be  one  of  an  arbitrarily  large  set  of  functions.  For  com¬ 
putational  efficiency,  we  assume  that  the  effect  of  the  material  and  strain  dependence  can  be  decomposed  by 
an  appropriate  choice  of  form  for  this  function.  This  can  be  done  by  selecting  a  form  expressed  in  terms  of 
a  set  of  coefficients  that  depend  solely  on  the  material  properties,  and  a  set  of  constant  basis  functions.  For 
simplicity,  the  following  linear  combination  case  is  selected: 

<j>(e,m)  =<)>(£,  c)  =  c^mJXjCe) +  ...+cm(m)xn(e)  =  ct  (m)  xf  (§)  •  02) 

In  this  form,  c  represents  the  vector  of  the  material  depended  coefficients  c,-,  and  X,  represents  the 
basis  functions  depending  only  on  strains  £  and  defined  at  a  total  of  n  distinct  points  distributed  over  the 
strain  space.  Equation  (12)  can  be  thought  as  being  an  interpolation  function  allowing  evaluation  of  (J>on 
points  other  than  the  ones  used  to  define  the  basis  functions. 

In  our  case,  125  (i.e.,  5  per  axis)  such  points  were  selected  to  be  uniformly  distributed  in  the  cartesian 
representation  of  the  components  of  the  strain  space.  Figure  14  shows  the  cuboidal  distribution  of  the  points 
where  the  basis  function  for  our  case  is  defined. 

To  have  a  unique  number  characterizing  each  one  of  the  points  where  the  strains  are  evaluated,  three 
indices  04),  and  q  have  been  chosen  along  the  E  ,  E  axes,  respectively,  to  count  from  1  to  5  as  the 
(opq)  frame  of  reference  indicates  in  Fig.  14.  This  resu'lted'm  using  a  counter  i  for  the  points  associated  with 
these  indices  in  the  following  relationship: 

i  =  o  +  5(p-  1)  +25  (q-  1)  for  ,  o, p,  q  =  1 . 5.  (13) 

The  increments  along  the  three  orthogonal  strain  component  directions  in  this  space  are  equidistanced 
and  the  step  is  the  quotient  of  the  maximum  value  for  any  of  the  strains  obtained  from  the  finite  element 
analysis  of  the  structure  under  consideration  divided  by  4  (number  of  spaces  between  the  discrete  points  of 
evaluation). 
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The  basis  functions  have  the  property  of  being  equal  to  the  value  of  the  function  at  the  evaluation  point 
(in  this  case,  1),  and  zero  at  any  other  evaluation  point  of  their  domain.  A  consequence  of  this  is  that  the 
value  of  the  function  (j>  according  to  Eq.  (12)  at  an  evaluation  point  is  equal  to  the  coefficient  of  the  basis 
function  for  that  point.  This,  it  turns  out,  is  a  very  convenient  feature  for  the  solution  scheme  described  later 
in  this  section. 

Figure  14  gives  a  graphic  definition  of  how  a  basis  function  behaves.  The  five  successive  planes  have 
been  shifted  out  of  the  strain  space  and  have  been  drawn  such  that  the  vertical  axis  represents  the  value  of 
the  function  x  •  The  function  has  been  evaluated  at  the  point  defined  by  (o  =  4,/>  =  3,  q  =  3),  or  (i  =  64),  and 
as  can  be  seen,  the  value  of  the  function  is  zero  at  all  nodal  points  of  the  strain  space  except  for  the  point 
where  it  evaluates  to  1.  There  it  remains  the  same  regardless  of  the  evaluation  point 

By  letting  £.  signify  the  strain  vector  for  each  point  i  for  all  123  labeled  points  in  the  strain  space,  the 
basis  functions  with  these  properties  can  be  expressed  as: 


X;(e:)  J  for  =  1,2, ...,  125.  (14) 

J  0,i=/ 

Introducing  this  relation  to  Eq.  (12),  for  one  of  these  points  i  we  obtain 

<t»  (§,.  m)  =  c,  (m)  X,  (E,)  +  ...  +  ct (m) x,  (§,)  +  ...  +  cm  (m) x„  (e.) 

(j >(£,,  m)  =  Cj  (m)0  +  ...  +ci(m)  1  +  ...  +cm(m)C  (15) 

which  after  evaluation  yields: 


<|>  (ef,  m)  =  cf(m).  (16) 

Thus,  Eq.  (12)  may  be  viewed  as  an  interpolation  function  that  linearly  determines  the  value  of  the  dissipated 
energy  density  function  at  any  point  in  the  strain  space  as  long  it  lies  in  the  cuboidal  region  spanned  by  the 
125  points. 

The  problem  of  determining  the  analytical  expression  for<|>  is  now  reduced  to  determining  the  coeffi¬ 
cients  q,  such  that  the  family  of  functions  represented  by  Eq.  (12)  is  restricted  to  one  which  represents  a 
“best  fit”  of  experimental  data. 

For  this  purpose,  we  assert  that  the  total  energy  <I>  dissipated  throughout  the  volume  of  the  structure 
under  consideration  ( s0 )  for  each  loading  condition  Iq  can  be  computed  as  the  volume  integral  of  the  dissi¬ 
pated  energy  function  over  the  spatial  extent  of  the  structure: 

=  j<|>(e,  c)dx.  (17) 

di0 

By  virtue  of  Eq.  (12),  this  relation  reduces  to 

=  f  (c,(m)x1(e(/0^)) +  -  +  cm(m)xn(e(/0.?)))4j  =  } ci(m)xl(§(/0,x))4  (18) 

ds„  ds0 
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Fig.  14  —  Strain  space  representation  and  characteristic  representation  of  an  indicative  basis  function 


Assuming  that  D  represents  an  experimentally  obtained  value  of  the  total  dissipated  energy  for  a  given  struc¬ 
ture  *0  and  a  given  loading  condition  Iq,  there  is  an  appropriate  selection  of  the  coefficients  c,-  such  that  the 
value  obtained  from  Eq.  (18)  approximates  this  measured  value.  This  can  be  expressed  by 

<l>(s0,l0) +e  =  D(s0,lQ),  (19) 

where  e  expresses  the  error  between  the  two  values  of  the  dissipated  energy. 
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3.4  Solution  and  Optimization  Scheme 

For  the  case  of  the  particular  specimen  geometry  sq  and  the  particular  loading  case  that  corresponds 
to  a  corresponding  experimental  point  p,  the  integral  of  Eq.  (18)  can  be  equivalenced  by  a  discretized  form 
indicating  the  sum  of  the  dissipated  energy  of  each  of  the  150  elements  of  the  specimen  structure.  Equation 
(19)  can  now  be  rewritten  as: 

110 

JJciXi(?e)V’e  +  eP  =  &.  (20) 

e  =  1 

where  Vg  represents  the  volume  of  the  element  e  for  the  loading  point  p. 

Since  this  represents  one  of  many  relations  we  can  write  for  different  increments  of  the  index  1  corre¬ 
sponding  to  different  loading  conditions,  the  set  of  these  relations  can  now  be  represented  in  a  matrix  form: 

[X]c  +  e  =  d,  (21) 

110 

where  [X] ,  c,  and  d  matrices  have  elements  described  by  X,p  =  ^  X,-  (§^e)  ^ e ,  c,-  and  d,  respective¬ 
ly.  We  selected  every  third  point  from  each  loading  path  (17  experimental  points  per  loading  path)  for  all  15 
loading  paths  for  a  total  of  255  experimental  points,  and  d  is  constructed  as  follows: 

/=  [Ufi/,  [tf’jl7.  ••  Ufft*-..,  [tf’isl7].  with  p  =  1,4,  ...,49,  (22) 

where 

Ufd  T  =  [Dlh  D\  D\,  . . .,  Z>49,]  .  (23) 

where  the  elements  are  obtained  from  Eq.  (10)  for  each  loading  path  /. 

Equation  (21)  represents  an  overdetermined  system  of  255  equations  with  125  unknowns.  However, 
the  requirement  that  the  dissipated  energy  is  a  monotonically  increasing  function  of  the  strains  imposes  a 
new  set  of  constraints  that  can  be  now  be  expressed  in  the  following  inequality: 

[M)cZ0.  (24) 

The  array  [M]  has  125  columns  and  100  rows,  since  it  represents  a  tuple  of  100  additional  constraints. 
It  is  defined  as: 


T 

where  Mc  for  c  =  1, ...,  100  are  row  matrices  corresponding  to  each  one  of  these  constraints. 

Thus  the  problem  of  determining  the  125  coefficients  c,-  is  reduced  to  the  solution  of  Eq.  (21)  under 
the  inequality  constraints  of  relation  (24).  This  combination  provides  a  total  of  355  relations,  and  it  repre- 


24 


MASTETAL. 


sents  a  classical  optimization  problem  with  inequality  constraints  where  the  objective  is  to  minimize  the 
function  (objective  function)  ||f||  such  that  both  Eqs.  (21)  and  (24)  are  satisfied.  This  is  a  standard  problem 
in  quadratic  programming  and  is  readily  solved  using  well-established  numerical  techniques  (see  for  exam¬ 
ple  Gill,  Murray,  and  Wright  (1981)).  Appendix  C  gives  a  vector  space  graphical  representation  of  the  enti¬ 
ties  involved  in  the  optimization  scheme. 

In  light  of  the  preceding  considerations,  the  dissipated  energy  is  now  fully  determined  at  all  125  points 
of  the  strain  space  since,  as  discussed  previously,  the  computed  coefficients  represent  the  values  of  the  dis¬ 
sipated  energy  at  those  points.  The  dissipated  energy  can,  therefore,  now  be  computed  at  any  point  in  any 
material  that  has  been  tested.  The  strains  determine  the  particular  subregion  in  the  cube  spanning  the  strain 
space,  such  that  a  linear  interpolation  performed  between  the  values  of  the  dissipated  energy  at  the  surround¬ 
ing  point  establishes  the  value  of  the  dissipated  energy  at  the  intermediate  point. 

As  a  consequence  of  this,  we  can  spatially  map  the  absorbed  energy  over  the  specimen  for  each  ex¬ 
perimental  point.  Figure  15  shows  the  evolution  of  the  dissipated  energy  distribution  for  the  experimental 
points  corresponding  to  20%,  40%,  and  60%  of  the  total  load  (points  10, 20,  and  30)  in  association  with  the 
total  dissipated  energy  evolution  curve  for  the  specimen  associated  with  loading  path  1 1 ,  repetition  1 .  These 
distributions  were  obtained  by  using  the  structural  simulator  described  in  Section  4.  In  these  distributions, 
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the  darker  the  shade  of  gray,  the  higher  the  value  of  the  dissipated  energy  density.  From  this  figure  it  can  be 
deduced  that  for  the  low  load  level  corresponding  to  20%  of  the  total,  the  total  dissipated  energy  evolution 
curves  shows  that  almost  the  entire  specimen  is  in  the  linear  region  because  of  the  very  small  value  of  the 
dissipated  energy,  and  the  corresponding  distribution  shows  indeed  that  only  a  very  small  region  on  the  up¬ 
per  right  comer  entered  into  the  nonlinear  domain  with  very  small  values.  As  the  load  increases,  the  total 
dissipated  energy  of  the  specimen  increases  and,  as  expected,  more  areas  in  the  specimen  enter  into  the  non¬ 
linear  region,  while  the  old  ones  intensify  their  absorption  levels. 

4.  APPLICATIONS  ON  NAVAL  STRUCTURES 

In  this  section  we  develop  the  idea  of  a  structural  response  simulator  and  show  how  the  computed 
energy  dissipation  functions  can  be  used  within  the  simulator  to  design  structural  components  and  predict 
various  aspects  of  their  mechanical  behavior. 

Because  dissipated  energy  is  taken  as  a  measure  of  internal  structural  damage,  one  use  of  a  simulator 
would  be  for  structural  designers  to  examine  simulated  spatial  dissipated  energy  distributions  to  form  selec¬ 
tion  criteria  for  evaluating  the  structural  response  of  structures  at  a  particular  state  of  the  design  process. 
Material  system  designers  could  also  use  a  simulator  to  generate  spatial  dissipated  energy  maps  from  which 
they  could  determine  the  material  response  of  a  composite  material  system,  and  then  select  sets  of  manufac¬ 
turing  parameters  to  optimize  the  damage  tolerance  of  new  systems.  Using  the  simulator  to  analyze  material 
softening  distribution  patterns  could  also  assist  in  selecting  a  material  system  for  a  particular  structural  de¬ 
sign.  Finally,  a  structural  response  simulator  could  even  be  used  to  qualify  material  systems  by  using  far 
more  realistic  qualification  criteria  than  are  currently  employed,  e.g.,  strength  or  fracture  toughness,  which 
only  indirectly  and  partially  assess  some  qualitative  aspects  of  the  materials.  Indeed,  it  is  shown  in  this  sec¬ 
tion  that  higher  strength  does  not  necessarily  imply  better  damage  tolerance. 

4.1  Structural  Response  Simulator 

To  provide  the  engineering  community  with  a  tool  for  addressing  the  types  of  problems  discussed 
above,  a  structural  response  simulator  was  created  that  takes  user-specified  geometries,  materials,  and  load¬ 
ing  conditions  as  inputs  and  creates  spatial  material  distributions  of  local  stiffness  loss  using  the  experimen¬ 
tally  obtained  dissipated  energy  functions.  Figure  16  shows  the  architecture  of  the  structural  response 
simulator  in  the  configuration  used  to  generate  the  results  presented  in  this  section.  It  has  since  evolved  into 
a  much  more  sophisticated  form  and  its  design  now  includes  such  advanced  functionality  as  user-simulator 
interaction  modelling  via  artificial  neural  nets  and  other  techniques. 

The  Fig.  16  block  diagram  shows  the  structure  of  the  simulator  in  terms  of  component  subsystems. 
The  user  interface  was  created  by  writing  a  library  of  function  calls  using  the  “PATRAN  Command  Lan¬ 
guage”  (PCL )  that  allows  the  user  to  specify  the  geometry,  the  material,  and  the  loading  for  the  structure  of 
interest.  Visualization  capabilities  are  provided  through  PDA’s  “PATRAN”  2.4  solid  modeling  package.  A 
finite  element  idealization  of  the  structure  is  achieved  through  appropriate  use  of  PCL  and  all  relevant  model 
information  is  stored  in  “PATRAN”  neutral  files  on  the  mass  storage  memory  resources  of  the  computer  sys¬ 
tem  hosting  the  simulator,  which  is  currently  a  Silicon  Graphics  Iris  4D/70GT  with  16  Mb  of  RAM  and  1.5 
GB  of  mass  storage  memory.  The  actual  finite  element  analysis  is  performed  by  using  the  “ABAQUS”  finite 
element  code  on  a  CRAY  X-MP-24  and  assuming  linear  elastic  behavior.  The  objective  of  the  finite  element 
analysis  is  to  generate  the  strain  field  in  the  structure  for  the  initial  elastic  state.  The  data  required  by 
“ABAQUS”  are  provided  by  the  “PATABA"  process  that  translates  the  “PATRAN”  neutral  files  to  an 
“ABAQUS”  input  file.  The  output  from  “ABAQUS”  is  stored  in  a  file  that  is  subsequently  converted  to  a 
strain  file  and  a  displacement  file  in  “PATRAN”  format  by  the  “ABAPAT”  process.  The  displacement  files 
are  identified  by  the  suffix  “dis”  and  the  strain  files  have  the  suffix  “.els.”  At  this  stage,  the  strains  and  dis- 
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Fig.  16  —  Block  diagram  of  the  architecture  of  the  structural  simulator  and  computational  resource  al'ocation 


placements  can  be  visualized  from  “PATRAN”  or  be  fed  into  the  dissipated  energy  module  “ABSORB” 
which  runs  on  a  Symbolics  3600  and  generates  a  spatial  energy  absorption  map  by  associating  a  value  of 
dissipated  energy  with  each  strain  combination  at  every  spatial  evaluation  point.  The  required  dissipated  en¬ 
ergy  values  are  computed  by  interpolation  using  the  coefficients  in  the  file  “rmt.x.y.pred”  that  were  produced 
by  the  process  depicted  in  Fig.  12.  Finally,  the  white  arrows  that  connect  the  user  interface  module  with  the 
other  simulator  modules  signify  user  control  over  these  modules  through  both  the  user  interface  and  the  net¬ 
working  TCP/IP  utilities  that  control  intermachine  data  transfer. 

42  Ship  Mast 

One  of  the  design  goals  is  to  investigate  the  merits  of  using  composite  materials  for  ship  masts.  The 
primary  reasons  composite  masts  are  being  considered  are  weight  reduction  and  favorable  radar  return  char¬ 
acteristics.  Traditional  designs  for  metal  masts  were  based  on  the  assumption  of  linear  elastic  behavior  and 
stress-based  damage  criteria;  namely,  that  internal  load-induced  damage  was  assumed  to  occur  whenever 
the  stress  exceeded  a  specified  maximum  stress  set  at  2.5  times  the  yield  stress  for  the  material  used.  As 
mentioned  previously,  such  a  simplistic  view  of  failure  behavior  is  not  justified  in  composites,  and  a  more 
realistic  damage  model  such  as  the  one  being  proposed  in  this  study  is  required.  The  enhanced  understand¬ 
ing  provided  by  such  a  model  should  lead  to  lower  safety  factors  and,  therefore,  more  efficient  mast  designs. 
The  structural  response  simulator  together  with  a  database  of  computed  dissipated  energy  functions  repre¬ 
sent  a  unique  medium  to  achieve  the  efficient  design  of  composite  masts. 

The  DDG-51  class  of  frigates  was  selected  as  a  candidate  ship  for  a  representative  mast  installation. 
The  geometry  of  the  mast  is  that  of  a  tripod  mast  with  a  rectangular  box  beam  and  two  struts  with  two  an¬ 
tenna  platforms  attached  to  it  and  the  associated  deckhouse  as  shown  in  Fig.  17. 

The  design  loads  selected  for  this  case  were  extracted  from  the  “DDS-170”  manual.  For  demonstra¬ 
tion  purposes,  a  wind  load  of  30  psf  (lb/ft2)  corresponding  to  a  90  kn  head  wind  was  selected.  In  addition, 
an  inertial  load  of  0.6  g’s  applied  downwards  to  the  center  of  gravity  of  the  mast,  or  0.2  g’s  applied  upwards 
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(a)  x  (b) 


Fig.  17  —  Rendered  view  of  (a)  the  mast  structure  and  (b)  wireframe  view  of  the  mast  structure  with  the 
view  plane  used  for  displaying  die  dissipated  energy  density  distributions 

on  the  center  of  gravity  of  the  mast.  A  500  lb  antenna  load  was  assumed  to  be  applied  at  the  upper  platform 
and  a  1000  lb  antenna  load  applied  at  the  lower  platform  of  the  mast  were  also  considered.  Only  the  case  of 
head-on  wind  load  is  presented  here,  however. 

The  materials  used  for  this  simulation  were  item  007  (Resin  3501-6,  Fiber  AS4,  Layup  +/-  60°)  and 
item  040  (Resin  PEEK,  Fiber  AS4,  Layup  +/-  60°)  as  described  in  Table  1.  These  materials  wane  selected 
to  establish  the  relative  merits  of  the  thermoset  3501-6  resin  material  vs  the  thermoplastic  PEEK  resin  ma¬ 
terial.  The  material  was  composed  of  130  successive  plies  of  0.0075  in.  each,  thus  accounting  for  a  total 
thickness  of  0.975  in. 

Figure  18  shows  the  dimensions  and  the  finite  element  mesh  for  the  simulation  of  the  mast  structure 
itself.  A  total  of  106  “ABAQUS”  “QUAD”  elements  were  used  for  the  finite  element  idealization.  In  addi¬ 
tion,  battle  damage  was  simulated  by  considering  a  5  inch  hole  at  the  right  side  of  the  foot  of  the  mast  as 
shown  in  Fig.  18. 

Figures  19  through  22  present  the  simulation  results  for  different  loading  and  material  choices  as  seen 
from  the  front  of  the  mast  in  the  legion  indicated  in  Fig.  17  with  a  rectangular  window.  All  figures  are  fringe 
p’ots  of  the  dissipated  energy  distribution  over  the  structure.  In  all  of  these  figures,  the  views  on  the  left  rep¬ 
resent  the  dissipated  energy  distribution  for  the  3501-6  thermoset  material,  and  the  views  on  the  right  rep¬ 
resent  the  dissipated  energy  distribution  for  the  PEEK  thermoplastic  material. 
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Fig.  18  —  Right  side  view  of  the  finite  element  meshing  and  dimensions  of  the  mast  (a);  front  view  of  the  mesh  (b);  cross 
section  of  the  mast  (c) 
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Fig.  19  —  Front  views  of  dissipated  energy  distributions  of  ship  mast  made  out  of  therm oset  material  007  ( AS4/ 
3501-6  [+/-600])  (left)  and  out  of  thermoplastic  material  040  (AS4/PEEK  [+/-600])  (right),  for  lx  of  the  total  load¬ 
ing  corresponding  to  front  wind 
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Fig.  20  —  Front  views  of  dissipated  energy  distributions  of  ship  mast  made  out  of  thermoset  material  007  (AS4/ 
3501-6  [+/-60p])  (left)  and  out  of  thermoplastic  material  040  (AS4/PEEK  [+/-600])  (right),  for  2x  of  the  total 
loading  corresponding  to  front  wind 
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Fig.  2 1  —  Front  views  of  dissipated  energy  distributions  of  ship  mast  made  out  of  thermoset  material  007  (AS4/ 
3501-6  [+/-60 ^1)  (left)  and  out  of  thermoplastic  material  040  (AS4/PEEK  (+/-60°J)  (right),  for  3*  of  the  total 
loading  corresponding  to  front  wind 
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Fig.  22  —  Front  views  of  dissipated  eneigy  distributions  of  ship  mast  made  out  of  thermoset  material  007  (AS4/ 
3501-6  [+A60°1)  0eft)  and  out  of  thermoplastic  material  040  (AS4/PEEK  [+/-600])  (right),  for  4*  of  the  total 
loading  corresponding  to  front  wind 
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The  effect  of  applying  an  increasingly  higher  load  corresponding  to  front  wind  is  manifested  as  an 
increasingly  higher  distribution  of  fringes  for  higher  levels  of  loading.  Each  one  of  these  figures  has  been 
plotted  for  lx,  2x,  3x,  and  4x  overload  levels  of  the  original  30  lb/ft2  loading  in  the  same  order  they  appear 
in  the  document. 

Another  characteristic  of  the  response  of  the  structure  due  to  the  increasing  load  is  the  effect  demon¬ 
strated  in  Fig.  23  where  the  maximum  dissipated  energy  is  plotted  as  a  function  of  the  load  magnitude.  That 
is,  after  a  certain  loading  level,  the  material  reaches  a  saturation  level,  where  the  dissipated  energy  in  the 
material  ceases  to  increase.  This  was  the  reason  the  distributions  of  overload  magnitudes  4  and  5  were  al¬ 
most  identical  and,  therefore,  it  seemed  appropriate  not  to  present  the  plotted  distributions  for  loading  mag¬ 
nitude  5.  Evidence  of  this  can  be  seen  by  comparing  the  dissipated  energy  distributions  for  loading 
magnitudes  3  and  4  in  Figs.  21  and  22. 

A  significant  observation  is  that  the  3501-6  material  seems  to  show  a  significantly  higher  propensity 
to  go  into  its  nonlinear  region  than  does  the  PEEK  material.  This  is  evident  from  the  higher  number  of  fring¬ 
es  and  their  higher  intensity  as  they  appear  on  the  3501-6  distributions  in  Figs.  19  through  22  compared  to 
those  of  the  PEEK  in  Figs.  19  through  22.  In  addition,  the  3501-6  material  also  has  the  tendency  to  dissipate 
higher  amounts  of  dissipated  energy  as  is  shown  by  the  saturation  levels  in  Fig.  23. 

As  is  apparent  from  all  these  figures,  the  mechanical  damage,  in  terms  of  dissipated  energy,  caused 
by  the  round  penetration  on  the  left  side  of  the  base  area  of  the  mast  is  negligible  compared  to  the  one  on 
the  middle  level  of  the  structure. 


Fig.  23  —  Distribution  of  the  maximum  dissipated  energy  of  the  ship  mast  structure  for  materials  007  (AS4/3501  -6  [+/-60p]) 
and  040  (AS4/PEEK  [+/  60°]) .  vs  the  loading  magnitude 
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43  In-Plane  Loader  Specimens 

The  specimen  used  in  the  IPL  tests  is  an  example  of  a  structural  object.  It  has  structural  loads,  i.e.,  the 
loads  applied  by  the  IPL  actuators.  It  has  a  structural  geometry,  i.e.,  the  geometry  of  the  specimen.  And  it 
has  a  structural  material,  i.e.,  the  material  under  investigation. 

At  any  point  within  the  structural  object,  strains  are  induced  that  act  upon  the  material.  In  this  instance, 
strains  induced  by  the  actuator  loads  act  upon  the  specimen  material.  The  strains  are  occasionally  the  same 
at  different  points  of  the  specimen,  but  in  general  they  are  not.  The  strains  act  upon  the  specimen  material, 
changing  it  on  occasion.  When  the  strains  do  change,  the  material  shows  the  change  by  exhibiting  a  different 
elastic  behavior  suggesting  reduction  in  local  stiffness  (softening). 

The  experimentally  determined  dissipated  energy  is  a  measure  of  material  softening.  It  is  a  nonlinear, 
monotonic  function  of  the  strains  that  act  upon  the  material.  As  such,  the  material  exhibits  behavioral  chang¬ 
es  that  may  be  unexpected. 

This  application  uses  the  IPL  structural  object  to  showcase  the  nature  of  the  unexpected  behavior. 

Continuing  in  the  spirit  of  Sections  2  and  3,  material  002  in  Table  1  was  used  as  the  straw  person  in 
this  application.  With  a  layup  construction  of  [(+30  /-30)4]s,  this  laminated  composite  material  has  plies 
constructed  from  AS1  carbon  fibers,  and  3501-6  epoxy  resin. 

Load  icons  were  employed  in  the  figures  of  this  application  to  facilitate  recognition  of  the  particular 
load  combination  being  applied  to  a  specimen.  Five  graphic  symbols  were  used  in  combination  to  produce 
these  icons.  They  are:  a  clockwise  arrow  to  indicate  clockwise  motion  of  the  movable  boundary  of  the  spec¬ 
imen;  a  counterclockwise  arrow  to  indicate  the  opposite  motion;  an  upwards-pointing  arrow  to  indicate  mo¬ 
tion  of  the  points  of  the  movable  boundary  parallel  to  the  y-axis;  a  downwards  arrow  for  the  opposite 
motion;  and  a  right-pointing  arrow  to  indicate  motion  of  the  points  of  the  movable  boundary  parallel  to  the 
x-axis. 


The  interpretation  of  a  particular  combination  of  the  graphic  symbols  is  that  the  corresponding  com¬ 
bination  of  IPL  applied  motion  is  applied  proportionally  to  the  specimen’s  movable  boundary.  The  speci¬ 
men,  caught  between  its  movable  and  immovable  boundaries,  deforms  with  a  compatible  motion.  Loads  that 
are  proportional  to  each  other  are  represented  by  the  same  icon.  The  particular  magnitude  of  the  proportional 
load  is  indicated  individually,  or  by  group  where  applicable. 

That  specimen  dissipated  energy  is  a  continuous  nonlinear  monotonic  Junction  of  load  displacement 
magnitude  is  an  experimentally  determined  fact.  The  material  dissipation  energy  is  forced  to  be  a  continuous 
nonlinear  monotonic  function  of  strain  magnitude.  How  this  was  done  is  explained  in  Section  3;  the  point 
being  that  though  the  specimen,  undeniably,  is  nonlinear  monotonic,  the  computed  material  dissipation 
function  behaves  that  way  by  analytic  design. 

Figure  24  d<  :s  the  loading  case  for  load  path  1 1  (lp- 1 1).  The  load  magnitudes  shown  are  20, 40, 
and  60  percent  of  the  maximum  boundary  motion  applied  during  the  experiment  for  this  particular  speci¬ 
men.  Figure  24  also  illustrates  the  monotonic  nature  of  dissipated  energy  with  respect  to  strains.  Dissipated 
energy  never  decreases  with  increasing  magnitudes  of  loading.  It  is  also  apparent  here  that  the  dissipated 
energy  may  not  increase  after  reaching  a  certain  magnitude  of  strain,  and  that  a  minimum  magnitude  of 
strain  is  required  before  the  dissipated  energy  can  become  greater  than  zero.  This  observed  behavior  corre¬ 
sponds  precisely  to  what  we  mean  by  the  term  a  nonlinear  monotonic  Junction  of  strain  magnitude. 
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Fig.  24  —  Distributions  of  dissipated  energy  on  the  deformed  specimen  for  loading  path  11,  material  002 
(AS 1/3 506-1  [+/-  30°])  for  three  different  magnitudes  of  loading 
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To  reiterate,  that  specimen  dissipated  energy  is  a  continuous  nonlinear  Junction  of  load  combination 
is  an  experimentally  determined  fact.  In  computing  the  material  dissipation  energy  function  from  the  spec¬ 
imen  dissipated  energy,  the  only  analytic  control  exercised  is  to  ensure  continuous  behavior  with  respect  to 
the  strain  combination.  Consequently,  the  nonlinear  effect  of  strain  combination  as  illustrated  in  Figs.  25, 
26,  and  27  is  closely  related  to  the  experimentally  observed  behavior. 

Figures  25, 26,  and  27  all  depict  the  same  set  of  the  15  different  load  combinations  used  in  obtaining 
the  experimental  results.  Figure  25  is  for  the  case  of  20%  of  full  load  magnitude,  while  Fig.  26  is  for  40%, 
and  Fig.  27  for  60%.  The  effect  of  opening  vs  closing  motion  (lp-14  vs  lp-2),  which  may  be  likened  to  ten¬ 
sion  vs  compression  behavior,  is  quite  distinctive.  At  the  20%  level,  both  cases  are  indistinguishable.  At  the 
40%  load  level,  the  opening  or  tension  case  has  developed  an  extensive  hot  area  from  the  notch  tip  out  along 
the  +/-  30°  fiber  directions  all  the  way  to  the  edge  of  the  specimen,  while  the  closing  or  compression  case 
has  changed  only  slightly.  At  the  60%  load  level,  the  situation  has  changed  dramatically  again.  The  com¬ 
pression  case  has  developed  a  hot  zone  of  the  same  configuration  as  that  for  tension.  But  while  the  zone  for 
opening  is  diffuse,  that  for  closing  is  concentrated  at  the  notch  and  at  the  edge. 

The  effect  of  opening,  sliding,  and  closing  load  combinations  (lp-14,  lp-8,  lp-2),  demonstrates  clearly 
the  distinction  between  linear  effect  and  nonlinear  effect  The  opening  and  sliding  combination  is  lp- 1 1 .  Lp- 
14  and  lp-8  both  show  extensive  zones  of  nonzero  dissipated  energy  in  the  lower  half  of  the  specimen.  A 
linear  effect  would  show  dissipated  energy  as  a  proportional  sum  of  the  dissipated  energy  from  cases  lp-14 
and  lp-8.  Clearly,  lp- 1 1  shows  no  nonzero  dissipated  energy  in  the  lower  half  of  the  specimen.  The  effect  is 
typical  of  nonlinearity.  The  sliding,  closing  combination,  lp-5,  exhibits  the  same  nonlinear  effect 

The  general  assertion  is  that  “if  the  IPL  experiments  were  conducted  using  another  specimen  shape, 
then  the  material  dissipation  energy  function  would  be  the  same.”  In  other  words,  the  material  behavior 
should  be  independent  of  the  structure  or  of  the  loads  that  may  be  applied  to  that  structure.  If  the  preceding 
discussion  had  been  about  this  alternative  experimental  specimen,  then  the  different  specimen  shape  would 
have  produced  different  strain  distributions  and  consequently  different  material  dissipated  energy  distribu¬ 
tions.  The  same  general  continuous,  nonlinear,  monotonic  nature  of  the  material  dissipation  energy  function 
would  have  been  noted,  but  their  distributions  within  the  specimen  would  have  been  different.  If  a  strong 
relationship  had  been  drawn  between  the  specimen,  its  loadings,  and  the  consequent  distributions  of  material 
dissipated  energy  from  the  cases  presented,  then  a  different  strong  relationship  would  have  been  drawn  from 
consideration  of  the  results  for  the  alternative  specimen. 

The  danger  is  in  attempting  to  infer  material  behavior  from  specimen  behavior,  without  regard  to  how 
the  strains  are  distributed  within  that  particular  specimen. 
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Fig.  25  —  Distributions  of  dissipated  energy  on  the  deformed  specimen  for  load  magnitude  2.  material  002 
(AS1/3506-1  [+/-  30°]),  for  each  of  the  load  paths  1-15 
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Fig.  27  —  Distributions  of  dissipated  energy  on  the  deformed  specimen  for  load  magnitude  4,  material  002 
(AS1/3506-1  [+/-  30°]),  for  each  of  the  load  paths  1-15 
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4.4  Submarine  Hull 

An  underwater  structure,  such  as  a  submarine,  is  subject  to  a  vast  range  of  loads.  These  loads  vary 
from  underwater  explosions,  underwater  wave  actions,  maneuvering  loads,  depth  loads,  to  a  whole  variety 
of  other  types  of  loads.  The  survivability  of  the  structure  under  these  conditions  is  of  {Mime  importance. 

Since  their  invention,  submarines  have  mainly  been  constructed  from  metals,  but  today,  composites 
are  being  considered  as  a  viable  alternative.  Dissipated  energy  offers  a  means  of  mapping  out  or  hot  spotting 
areas  of  concern  on  composite  submarine  structures.  The  structural  loads,  structural  geometry,  and  the  struc¬ 
tural  material  act  to  produce  strains  throughout  the  structure  that  may  induce  changes  in  the  material.  Dis¬ 
sipated  energy  as  a  function  of  strain  determines  when  a  material  will  change  from  one  elastic  st  'e  to 
another,  and  is  an  important  means  of  measuring  both  the  extent  of  the  material  changes  and  the  structural 
survivability. 

In  this  application,  we  consider  a  generic  submarine  and,  in  particular,  analyze  the  cylindrical  hull  sec¬ 
tion.  The  actual  dimensions  used  are  those  of  a  generic  dry  dock  shelter.  The  cylindrical  shell  idealization 
was  selected  as  a  simple  showcase  for  dissipated  energy  mapping. 

The  finite  element  model  used  is  shown  in  Fig.  28  and  consists  of  256  shell  “QUAD”  elements.  The 
cylinder  that  was  analyzed  was  100  inches  long,  80  inches  diameter,  and  1  inch  thick.  The  boundary  loading 
conditions  were  those  induced  by  rigid  body  plane  rotations  of  the  ends  of  the  cylinder.  These  loading  con¬ 
ditions  where  chosen  as  a  nominal  representation  of  the  kinds  of  loads  that  might  be  carried  through  to  the 
cylinder  by  the  rest  of  the  submarine  in  response  to  bending  loads  about  the  vertical  or  y-axis.The  shell  ma¬ 
terial  was  a  laminated  graphite  thermoset  epoxy.  The  layers  were  oriented  at  +/-  60  °  to  the  longitudinal,  or 
z-axis.  The  fibers  were  AS4  and  the  resin  was  a  3501-6  type  epoxy. 

The  finite  element  code  “ABAQUS”  was  used  to  determine  the  interior  strain  field.The  strains  as  re¬ 
ported  back  by  the  code  are  for  the  Gauss  integration  points  used  by  “ABAQUS”  to  determine  the  shell  prop¬ 
erties.  Three  integration  points  were  used  in  this  analysis,  one  point  at  the  center  of  the  laminate  and  two 
points  at  about  an  eighth  of  the  thickness  in  from  each  side.  In  this  case,  the  dissipated  energy  is  not  com¬ 
puted  on  a  per  ply  basis,  but  rather  on  a  smeared-out  point  basis,  or  a  partial  dissipated  energy  number  at 
an  integration  point.  Gaussian  integration  is  then  used  to  compute  a  total  dissipated  energy  through  the 
thickness.  Both  types  of  dissipated  energy  are  illustrated  in  the  following  figures  where  the  partial  dissipated 
energy  at  the  outer  part  of  the  shell  is  shown  in  Fig.  28  and  total  dissipated  energy  in  Fig.  29. 

Figure  29  depicts  the  dissipated  energy  maps  for  overloads  of  7X,  9X,  1 1 X,  1 3X  and  1 5X.  The  values 
shown  are  for  partial  dissipated  energy  close  to  the  outer  side  of  the  shell.  As  can  be  seen  in  Fig.  29,  the 
boundary  of  the  hot-spotted  zone  (higher  values  of  dissipated  energy)  at  7X  overload  is  practically  the  same 
as  at  15X  overload.  An  observer  may  interpret  this  in  several  ways;  however,  it  is  reasonable  to  surmise  that 
the  damage  as  indicated  by  the  hot-spotting  is  not  going  to  grow  unstably,  and  that  regardless  of  the  point 
of  view  taken  as  to  the  meaning  of  the  damage,  the  material  degrades  only  so  far  and  thereafter  gets  no 
worse. 

Figure  30  depicts  the  distribution  of  dissipated  energy  from  below  the  inside  surface,  the  outside  sur¬ 
face,  and  throughout  the  thickness,  i.e,.  partial  vs  total  viewpoints.  The  partial  viewpoints  show  both  simi¬ 
larities  and  differences.  In  both,  the  extent  of  the  hot  zone  is  roughly  the  same;  however,  each  exhibits 
markedly  different  distributions  within  the  zone. 
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Fig.  29  —  Dissipated  eneigy  distributions  at  the  outer  surface  of  the  shell  structure,  for  loading  magnitudes  of  7,9,1 1,13, 
and  15  times  overload 


The  total  dissipated  energy  viewpoint  shows  a  larger  hot  zone  than  either  of  the  partial  dissipated  en¬ 
ergy  maps.  It  also  shows  a  rather  constant  value  of  dissipated  energy  over  the  area  of  the  hot  zone. 

From  this  an  important  inference  may  be  drawn,  that  the  nonlinear  effect  of  strain  on  material  loss  of 
local  stiffness,  or  change,  may  not  be  characterizable  strictly  on  the  basis  of  a  partial  or  a  total  viewpoint. 

5.  A  GENERAL  PROCEDURE  FOR  DETERMINING  THE  DISSIPATED  ENERGY  DENSITY 
FUNCTION  AND  A  DESCRIPTION  OF  MATERIAL  CONSTITUTIVE  BEHAVIOR 

The  approximate  procedure  for  determining  the  dissipated  energy  function  <j>  (e)  discussed  in  Sec¬ 
tion  3  is  currently  included  in  the  response  simulator.  As  presently  formulated,  this  procedure  has  the  fol¬ 
lowing  shortcomings:  first,  the  required  strain  fields  are  computed  by  assuming  linear  elastic  behavior  and 
by  using  the  elastic  constants  corresponding  to  the  "virgin”  material,  rather  than  by  using  the  proper  nonlin¬ 
ear  constitutive  behavior,  and  second,  the  procedure  is  formulated  only  in  terms  of  a  particular  piecewise 
linear  collocation  set  of  basis  functions  for  representing  <(>(£),  rather  than  in  terms  of  arbitrary  basis  func¬ 
tions. 
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Fig.  30  —  Dissipated  energy  distribution  for  (a)  the  inside  surface  of  the  shell,  (b)  the  outside  surface  of  the  shell  and 
(c)  integrated  through  the  thickness. 

Both  of  these  issues  are  addressed  in  the  present  section  and  the  underlying  assumptions  involved  in  the  ma¬ 
terial  model  are  carefully  laid  out. 

To  determine  the  dissipated  energy  function,  we  regard  the  composite  as  being  composed  of  either  a 
single  mechanically  equivalent  homogeneous  anisotropic  material,  or  a  collection  of  layers  of  varying  ori¬ 
entations  of  such  materials.  Provided  that  the  applied  loads  are  either  quasi-static  or  dominated  by  low  tem¬ 
poral  frequencies,  these  homogenization  procedures  should  provide  acceptable  models  since  it  is  expected 
that  the  wavelengths  corresponding  to  the  spatial  variation  of  the  stresses,  strains,  etc.  should  be  large  com¬ 
pared  to  the  microstructural  characteristic  lengths.  If  this  is  not  the  case,  events  on  the  micro  scale  can  take 
on  a  predominant  role  and  cause  significant  inaccuracies. 

We  also  postulate  that  the  material  constitutive  behavior  can  be  completely  described  by  specifying 
stress  as  a  function  of  only  strain  and  a  set  of  internal  state  variables,  i.e., 

g  =  Cft;,  e),  (26) 
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where  C(4*§)  is  a  nonlinear  fourth  order  tensor  function  and  £  is  the  state  variable  vector.  The  form  of  Eq. 
(26)  thus  obviously  precludes  dealing  with  materials  that  exhibit  marked  strain-rate-dependent  behavior.  Fi¬ 
nally,  we  ignore  any  residual  stresses  or  strains  which  may  exist  when  the  loads  are  removed  such  as  those 
induced  by  the  curing  process. 

The  form  of  th:  above  constitutive  relation  implies  a  restricted  implicit  type  of  path  dependency  in 
the  sense  that  any  path  dependent  behavior  results  only  from  internal  state  transitions  that  are  reflected  by 
changes  in  the  state  variables  and  not  from  an  explicitly  stated  load  history  dependence.  This  section  con¬ 
siders  this  issue  in  more  detail.  It  suffices  to  state  here  that  for  the  most  part  the  material  behaves  in  a  manner 
similar  to  that  of  a  hyperelastic  material  (a  nonlinear  elastic  material  with  a  strain  energy  potential)  as  long 
as  the  state  variables  £  remain  unchanged.  Thus,  for  example,  when  a  material  point  is  loaded  up  to  a  spec¬ 
ified  strain,  the  subsequent  structural  response  is  dependent  solely  on  the  current  strains  and  is  independent 
of  deformation  history  provided  no  unloading  takes  place.  Upon  unloading,  the  material  behaves  elastically 
in  that  no  further  internal  material  damage  takes  place,  however,  a  state  transition  occurs  when  the  unloading 
commences,  and  it  is  this  state  transition  that  is  responsible  for  any  memory  the  material  might  have  of  its 
previous  loading  history.  As  mentioned  in  the  introduction,  this  material  model  and  its  particular  type  of  path 
independence  is  similar  to  the  model  employed  by  Schapery  (1987, 89a,  9b,  90).  Moreover,  it  is  consistent 
with  the  observations  of  Lambom  and  Schapery  (1988)  who  observed  a  type  of  path  independence  for  a 
range  of  tensile  and  torsional  loads  in  several  graphite-epoxy  systems. 

The  remainder  of  this  section  describes  a  general  procedure  for  determining  the  dissipation  density 
function  using  boundary  force  and  displacement  data  obtained  from  IPL  tests.  The  procedure  is  a  deconvo¬ 
lution  process  in  the  sense  that  the  observations  reflect  both  the  effects  of  material  behavior  and  specimen 
geometry,  and  that  the  geometric  effects  must  be  factored  out  to  obtain  information  relating  to  the  material 
alone.  We  also  present  a  derivation  of  the  precise  analytical  form  of  the  constitutive  equations. 

5.1  Caveats 

In  the  course  of  our  analysis  we  have  introduced  a  number  of  assumptions  and  simplifications,  some 
having  already  been  mentioned.  Nevertheless,  we  itemize  them  below  in  their  entirety  for  the  sake  of  com¬ 
pleteness: 

•  The  material  can  be  regarded  as  a  mechanically  equivalent  homogeneous  anisotropic  material  as  dis¬ 
cussed  above.  In  an  angle-ply  composite,  this  assumption  can  be  applied  either  on  a  ply-by-ply  basis 
or  to  the  collection,  as  appropriate. 

•  Loading  is  either  static  or  slowly  varying  in  accordance  with  the  considerations  already  discussed. 

•  The  material  behavior  can  be  represented  as  !?  =  C(^»e)  discussed  previously.  This  assumption  by 
definition  precludes  consideration  of  materials  that  exhibit  marked  viscous,  rate,  and  explicit  load  his¬ 
tory  dependencies.  Moreover,  given  this  assumption,  a  so-called  work  potential  (energy  per  unit  vol¬ 
ume),  v  (e) ,  can  be  defined  such  that  a  =  grade\y  (£,  e)  (Schapery  1987). 

•  The  total  energy  absorbed  by  the  material  during  loading  can  be  regarded  as  being  composed  of  the 
sum  of  a  reversible  (recoverable)  and  an  irreversible  (dissipative)  part.  The  reversible  component  is 
the  energy  that  would  be  recovered  if  the  material  were  to  unload,  whereas  the  irreversible  part  repre¬ 
sents  the  energy  which  is  dissipated  by  the  internal  damage  mechanisms.  The  latter  can  be  described 
by  a  dissipation  density  function  <j>  (§)  (dissipation  energy  per  unit  volume). 

•  The  constitutive  relation  is  continuous  in  both  stress  and  strain. 

•  The  deformations  are  sufficiently  small  so  that  the  infinitesimal  stress  and  strain  tensors  may  be  em¬ 
ployed. 
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•  Only  shell-like  structures  are  considered  so  that  the  stress  and  strain  fields  in  either  the  entire  material 
or  in  each  ply  are  two-dimensional  with  no  transverse  components,  i.e.,  the  stresses  and  strains  can  be 
regarded  as  vectors  having  the  form  g  =  ( C  ,  G„,  CT^)7 ,  and  E  =  (e^,  E^,  e  ,  where 
(T|,  £,  £)  is  a  coordinate  system  embedded  in  tne  malerial  or  layer  with  one  axis  along  the  fiber  direc¬ 
tion  as  shown  in  Fig.  31. 

•  Displacement  continuity  is  maintained  between  layers. 

It  should  be  stressed  here  that  no  explicit  assumptions  are  made  concerning  the  detailed  nature  of  the 
various  failure  processes.  Rather,  the  approach  is  to  introduce  a  minimal  number  of  hypotheses  that  we  feel 
are  in  accordance  with  physical  fact  and  readily  defensible.  Moreover,  many  of  the  above  restrictions  can  be 
relaxed  when  the  situation  warrants. 

SJ.  Analysis  Overview 

The  primary  issue  here  is  the  estimation  of  the  dissipated  energy  density  function  $  (§)  and  the  sub¬ 
sequent  computation  of  constitutive  behavior.  We  summarize  the  steps  below  and  provide  more  detailed  ex¬ 
planations  in  the  remaining  subsections. 

•  A  representation  is  chosen  for  the  dissipation  function  in  terms  of  a  set  of  m  basis  functions 

X  =  (xv  x2.  •••.  Xm)  and  an  initially  unknown  parameter  vector  c=  (clt  c2, ...» cm)  . 

•  A  uniform  set  of  loading  paths  in  displacement  space  is  selected  (15  for  each  material  or  15  for  each 
of  a  set  of  layup  configurations  for  angle-ply  composites),  boundary  forces  and  displacements  ((,  u) 
are  measured  at  50  equally  spaced  points  on  each  loading  path  as  described  in  Section  2,  and  the  en¬ 
ergy  imparted  to  the  specimen,  W ,  is  computed  for  each  observation  point  p.  If  the  composite  is  re¬ 
garded  as  a  single  equivalent  homogeneous  material  (in  the  case  of  an  angle-ply  composite,  this 
implies  that  each  layup  configuration  is  to  be  counted  as  a  separate  material),  then  this  procedure 
yields  n=750  values  of  W*7 per  material.  If  the  material  is  an  angle-ply  composite  and  we  wish  to  an- 


Fig.  31  —  Coordinate  system  embedded  in  the  shell  structure 
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alyze  the  constitutive  behavior  on  a  ply-by-ply  basis,  then  the  procedure  yields  n  =  750  X  np  values 
of  for  each  material,  where  np  denotes  the  number  of  layup  configurations. 

•  The  irreversible  portion  of  the  imparted  energy,  DP,  is  computed  for  each  observation  point  p  using 
the  relationship 


Tf  =  W*  -  ^  {f  •  /)  (27) 

where  i/are  the  boundary  displacements  and  tpare  the  reaction  forces. 

•  The  assertion  that  /^equals  the  integral  of  <}>  (g  (x)  )  over  the  volume  of  the  specimen,  where 

x  =  (t|,  £,  \)  is  applied  at  each  of  the  n  observation  points  p  and,  in  conjunction  with  a  represen¬ 
tation  of  the  constitutive  relation  in  terms  of  <j>  (e)  to  be  described  later,  the  strain-displacement  and 
equilibrium  relations,  and  any  appropriate  additional  constraints,  results  in  a  highly  overdetermined 
set  of  nonlinear  equations  for  the  m  components  of  the  parameter  vector  c  (a  system  of  n  x  m  equa¬ 
tions  with  n»  m). 

•  These  equations  are  solved  numerically  for  c  and  the  dissipation  density  function  <J)  (£)  computed. 

•  The  aforementioned  representation  of  the  constitutive  relation  is  used  in  conjunction  with  the  com¬ 
puted  values  of  c  to  determine  the  constitutive  behavior  o  =  C  (£,  §) . 


5.3  Representation  of  the  Dissipation  Function 

The  dissipation  density  function  0  (§)  is  represented  by  the  linear  combination 

<K§)  *  c-x(e)  (28) 

where  x  ( §)  is  a  vector  of  suitable  C°  (or  smoother  basis)  functions  over  the  3-space  (  ETin’  eCC’  *  and 

c  the  parameter  vector  to  be  determined.  This  representation  is  quite  general  in  the  sense  that  X(-  may  be 
chosen  more  or  less  arbitrarily,  for  instance,  as  interpolation  basis  functions  over  a  suitable  mesh  defined  on 
(e^,  Ej-j.,  E^)  in  which  case  any  of  the  usual  three-dimensional  finite  element  shape  functions  may  be 
used;  as  B-splines  or  cardinal  splines;  or  as  locally  defined  or  global  orthogonal  polynomials.  The  only  sig¬ 
nificant  restriction  on  the  form  of  the  basis  functions  is  that  positivity  of  the  quantity  <j>  (e)  must  be  assured 
to  maintain  agreement  with  physical  observation.  When  X  is  30  interpolation  basis,  the  components  of  the 
parameter  vector  c  are  simply  the  values  of  $  at  the  node  points;  however,  this  is  not  true  when  the  later 
two  representations  for  X  am  employed.  In  practice,  we  almost  always  chose  to  use  a  locally  defined  linear 
interpolation  basis  and  this  is  indeed  the  situation  described  in  Section  3. 

5.4  The  Constitutive  Relations 

In  this  section  we  present  an  overview  of  the  material  model  and  describe  in  detail  a  procedure  for 
computing  the  constitutive  behavior  in  terms  of  the  dissipated  energy  density  function  <{>  (g) . 

The  basic  premise  is  that  the  material  can  be  regarded  as  always  being  in  one  of  two  distinct  domains: 
namely,  either  purely  elastic  or  inelastic.  When  the  material  is  in  an  elastic  domain,  no  internal  damage  is 
presumed  to  occur;  this  results  in  behavior  that  is  load  history  independent  and  reversible  in  the  sense  that 
no  energy  dissipation  takes  place.  The  material  enters  an  inelastic  domain  during  loading  after  a  certain 
strain-dependent  threshold  is  crossed.  Here  progressive  internal  damage  leads  to  behavior  that  is  irreversible 
in  nature  with  energy  being  dissipated  by  the  various  internal  damage  mechanisms.  A  restricted  type  of  load 
path  independence,  whose  precise  nature  will  be  discussed  shortly,  is  applicable  in  this  domain. 
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One  way  to  model  this  type  of  behavior  is  to  postulate  the  existence  of  a  surface  T  (^)  such  as 
F  (g,  4)  =  0  in  strain  space  with  the  following  property:  namely,  that  whenever  the  strain  g  associated 
with  some  material  point  P  is  either  inside  the  surface  T  (^)  ,  or  outside  T  (£,)  with  decreasing  dissipated 
energy  (  d§  (g)  / d\e\  <  0),  then  P  is  in  an  elastic  domain'  When  g  is  outside  of  T  (4)  with  increasing 
dissipated  energy  (  flf0  (g) /d\e\  >  0),  P  is  in  an  inelastic  domain.  Here  the  vector  ^  Ts  a  vector  of  internal 
state  variables  that  can  be  thought  of  as  damage  parameters  and  the  surface  T  (^)'acts  as  a  threshold  that 
defines  the  elastic/inelastic  transition.  The  surface  T  (£)  thus  plays  a  role  that  is  loosely  analogous  to  that 
of  the  yield  surface  (in  stress  space)  in  classical  plasticity  theory  with  the  state  variables  £  in  the  role  of 
strain  hardening  parameters. 

The  model  we  used  is  based  on  the  above  considerations  with  the  additional  assumptions  that  there  is 
only  a  single  scalar  state  variable  t  that  takes  on  di  Crete  values  ^ ,  (j=0,l ...),  and  that  the  function 
F(e,%)  in  the  definition  of  T  (q)  be  identified  with  level  surfaces  of  the  dissipation  function  0 (g) .  A 
state  transition  takes  place  each  time  the  material  switches  from  an  inelastic  to  an  elastic  domain.  The  state 
variable  is  defined  by 


=  0 

for  material  in  the  initial  state,  and  by 


(29) 


^  =  max  0(g)  ,j  =  1,2,...  (30) 

otherwise,  where  max  0  (g)  represents  the  maximum  value  of  energy  dissipation  encountered  in  the  previ¬ 
ous  inelastic  domain.  Also,  the  function  F(§,^)  is  taken  in  the  form 

F( E,V)  =  ♦  (§)  (31) 

We  assert  that  the  state  variable  £  is  monotonically  increasing,  which  implies  that  the  threshold  sur¬ 
face  r  (4)  always  expands  in  strain  space  in  the  sense  that  the  elastic  domain  interior  to  the  surface 
T(^)  must  always  be  included  in  the  elastic  domain  interior  to  T  ( £/  )  .  Also,  the  parameters  that  occur 
in  the  description  of  the  constitutive  behavior  in  both  the  elastic  anJ  inelastic  domains  are  functions  of  the 
state  variable  ^ ,  thereby  implying  the  existence  of  a  denumerable  set  of  elastic  and  inelastic  domains  pa¬ 
rametrized  by  £/  alone.  This  leads  to  the  important  conclusion  that  the  only  memory  the  material  has  of  its 
past  is  via  the  state  variable  Q ,  which  implies  a  one  parameter  type  of  history  dependence. 


Some  aspects  of  the  inelastic  behavior  implied  by  this  model  are  illustrated  in  Fig.  32,  which  shows 
how  the  dissipation  density  function  and  an  arbitrary  stress  component  <5w  vary  with  g  alonp  three  arbitrary 
radial  paths  in  strain  space  in  the  two  inelastic  domains  bounded  by  the  surfaces  r° and  T  .respectively. 
Here  TJ  denotes  T  (q)  and  increasing  values  of  j  imply  increasing  values  of  the  associated  state  variable 


Inelastic  Behavior 

We  examine  the  inelastic  behavior  regime  first.  Let  us  assume  that  the  work  potential  y  (g)  can  be 
taken  as 


¥(§)  =  <&(§)  +0(§) 


(32) 
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Fig.  32  —  Dissipated  energy  density  as  a  function  of  (a)  strains,  and  (b)  corresponding  stress-strain  space 


where 


<*>(§)  =  ~  (a(e)  e) 

1  (33) 

We  chose  this  particular  form  for  Eq.  (33)  because  the  term  d>  (e)  is  equal  to  the  energy  density  re¬ 
covered  if  the  material  were  to  unload  linearly  from  its  present  state,  and,  as  is  shown  in  this  section,  is  en¬ 
ergetically  consistent  with  the  postulated  elastic  unloading  behavior.  It  is  obviously  not  the  only  choice  we 
could  have  made,  but  it  is  by  far  the  simplest. 

Since  g  =  gradEy  (Q,  e) ,  by  definition,  taking  the  gradient  of  Eq.  (32)  gives 

g  =  grad§<D  (e)  +  gradE(p  (e) 

W  v 

Equation  (34)  is  a  system  of  first  order  linear  partial  differential  equations  in  g  (e)  and  can  be  dealt 
with  routinely  the  method  of  characteristics  (Courant  and  Hilbert  (1962))  when  the  proper  initial  conditions 
are  specified  on  the  surface  T  ;  namely, 


P(§r)  =  [£*(8^,  (Ctj,  cx2))] 
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Here  (a,,  a2)  316  surface  coordinates  on  -  =  £  ai>  a2^  denotes  the  set  of 

e  €  r  (Q)  and  are  regarded  as  specified,  and  |_C*  (e  ,  (otj,  a2) )  J  is  the  known  matrix  representation 
of  the  constitutive  relation  for  the  previous  elastic  domain. 

In  general  we  cannot  find  explicit  closed-form  expressions  for  o  (§)  that  satisfy  the  initial  conditions 
of  Eq.  (35).  There  is,  however,  one  such  case  where  such  a  solution  can  indeed  be  found;  namely  the  case 
where  the  constitutive  relation  for  the  prior  elastic  domain,  i.e.,  Eq.  (35)  is  linear  in  E,  i.e., 

g(e)  =  [c* (e^,  (a,, <x2))]e,  (36) 

and  where  the  matrix  [c*  (E1*,  (<Xj,  a2) )  ]  is  constant  over  the  transition  surface  V  (tj)  ,  i.e., 

[Cf  (e1*,  (aj  ,  a2)  )  J  =  [Cf  (%f)  ]  on  P.  These  criteria  are  certainly  met  when  the  elastic  domain  cor¬ 
responds  to  that  of  the  “virgin”  material  assuming  it  to  be  linear  and  homogeneous  in  a  mechanical  sense. 
Thus,  taking  the  scalar  product  of  Eq.  (34)  with  £  and  using  the  definition  of  <I>  (e)  decouples  the  system 
of  Eq.  (34)  and  results  in  the  following  single  first  order  linear  partial  differential  equation  for  0(e); 


e  •  grad§0  -  20  =  - (e  •  grad£<p  (£)  ) ,  (37) 

with  initial  conditions 

o  (e1*)  =  ~  (§r  (^',  otp  a2)  [C*(?,  ctj,  a2)  ]  e1*  (%,  a,,  a2) ) .  (38) 

Equation  (37)  with  initial  conditions  of  Eq.  (38)  again  can  be  solved  by  means  of  characteristics  to 
yield  the  following  expression  for  O  (§) 


0(#favarp)  =  p2 


(O^)) 


f 1  d^jr) 

-{?**)• 


(39) 


evaluated  on  the  characteristic  lines.  The  above  integral  is  taken  along  a  characteristic  line,  that  are  defined 
by 


and  where  by  -r=  we  mean 


§  =  §  (^,  dj,  a2,  P)  =  Er  (^',  a,,  a2)  p. 


(40) 


^  (%,ava2)).  (41) 

Equation  (40)  can  be  regarded  as  a  coordinate  transformation  relating  the  coordinates  (d  j,  a2,  P) 
to  (e!>  ^  ■  To  obtain  CT(£)  we  also  need  the  inverse  transformation  (Ej,£2,  £3)  — >  (dj,  a2,  P) 

and  we  must  be  able  to  evaluate  the  gradients  grad£gr  (£)  ,  grad£P  (e)  ,  and  grad£Ot  (§)  ,  with 
a  =  (ar  a2) .  A  closed-form  expression  for  the  inverse  can  be  found  when  the  =  1,2)  are  taken 

as  the  cosines  of  the  angles  that  the  position  vector  |e  |  makes  with  the  E^  axes,  i.e., 
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a..  = 


y_ 

si 


Then  using  Eq.  (40),  Eq.  (42)  becomes 


°Sr  ~  |e|  • 


and  from  Eq.  (40),  we  get  for  P : 


(42) 


(43) 


(44) 


where  the  term  £  in  Eq.  (44)  is  now  understood  to  be  a  function  of  £  whose  values  are  obtained  by  substi¬ 
tuting  Eq.  (43)  into  e”  (§*,  a  j,  a2)  .  Equations  (43)  and  (44)  therefore  define  the  sought  after  coordinate 
transformation  (e^  e2,  e3)  (olf  a2,  p)  . 


j-rf 

The  gradients  [gradea  (e)  ] ,  grad t.e  (e) ,  and  gntdef3  (e)  are  now  readily  evaluated  by  using 
Eqs.  (40)  and  (42)  through '(44).  Thus, 


where 


grad?txy(e) 


1  0 
0  1 
0  0 


grad§er'(e)]  =  [grad(Xer'(e)]  [grad?a(e)] , 


(45) 


(46) 

(47) 


where  [grad aerV(e)]  = 
determined,  and 


'-“’JUk, 


is  known  since  the  transition  surface  is  already 


grad§p2(e)  =  ^-^(I-pfgradgE^e)]  )je,  (48) 

where  I  is  the  identity  matrix. 

The  stress  fieldo(e)can  now  be  obtained  by  substituting  the  expression  forO(4,  Ctj,  a2,  p)  inEq. 

(39)  into  Eq.  (34)  and  using  the  relationships  between  £  and  (^,  a(,  a2,  P)  as  given  in  Eq.  (40).  Thus 
a  (g)  is  given  by 
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(  0  r  p  - 
9  =  I?)  +  -=~2  (I -p [graded)]  )f § 

l  |e1  “  «J **  J 

“P2ji  [graded  (§)  ]grada 


Elastic  behavior 

In  the  elastic  regime  we  assume  that  the  material  can  be  modeled  as  a  particular  type  of  nonlinear  hy¬ 
perelastic  material;  namely  one  that  allows  for  no  initial  stresses  and  has  linear  stress-strain  behavior  along 
radial  paths  in  strain  space.  This  later  assumption  was  invoked  to  maintain  energetic  consistency  with  the 
representation  for  \|f  in  Eq.  (32)  so  that  the  eneigy  available  for  elastic  unloading  is  always  given  by  O  (e)  . 
This  assumption  also  agrees  with  repeated  observations  of  linear  unloading  behavior  during  the  IPL  tests. 

The  type  of  behavior  we  are  attempting  to  model  is  illustrated  in  Fig.  33  where  the  linear  behavior 
along  radial  paths  is  shown  for  two  elastic  domains  with  associated  transition  surfaces  I~° and  T1  (Fig.  33 
(a)  and  (b),  respectively).  The  important  point  to  note  here  is  that  the  elastic  constitutive  behavior  is  gener¬ 
ally  nonlinear;  it  is  only  guaranteed  to  be  linear  along  the  radial  paths  mentioned  above.  Also  note  that  the 


Fig.  33  —  Stress-strain  behavior  for  state  I"°  ( a )  and  state  F1  ( b ) 
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effective  elastic  stiffness  is  expected  to  decrease  with  increasing  values  of  the  state  variable  .which  in  turn 
correspond  to  increasing  values  of  energy  dissipation  as  illustrated  in  the  Fig.  33. 

Consistent  with  the  above  assumptions,  we  propose  a  constitutive  relation  in  the  elastic  domain  having 
the  form 


9  =  (C'(4»a1,a2))|E|,  (50) 

where  Ce  is  now  taken  to  be  a  vector  rather  than  a  matrix.  Equation  (50)  has  a  deceptively  simple  form; 
however  it  is  really  quite  general  and  easily  subsumes  the  usual  linear  form 

g  =  [A]e.  (51) 

That  this  is  indeed  the  case  is  readily  demonstrated  since  £  can  always  be  written  as 
£  =  e  (ctj,  a2)  |e|  along  some  path  with  £  being  the  associated  unit  vector.  Therefore  Eq.  (51)  assumes 
the  form  of  Eq.  (50)  with  Ce  given  as 

Ce(a,,a2)  =  [A]l(aj,  a2).  (52) 

For  the  initial  state  the  constitutive  behavior  is  assumed  to  be  given  by  the  equivalent  form  of  Eq. 
(51)  with  the  matrix  A  computed  using  the  “virgin”  material  properties.  For  subsequent  states 

i  =  1,2, ... ,  Eq.  (50)  is  used  where  stress  continuity  across  the  surface  I*  leads  to  the  following  ex¬ 
pression  for  the  vector  Ce  E) 


C' (%'*',{) 


lei" 


where  a  ( E  )  is  computed  by  evaluating  the  stresses  given  by  Eq.  (44)  on  the  surface  F. 


(53) 


We  stress  that  Eq.  (51)  with  Ce  given  by  Eq.  (53)  usually  results  in  elastic  constitutive  behavior  that 
is  inherently  nonlinear  because  Ces  (£,  §)  is  generally  not  of  the  form  given  in  Eq.  (52). 


5.5  Equations  for  the  Determination  of  the  Parameter  Vector 

The  irreversible  portion  R1  of  the  total  energy  imparted  to  the  specimen  via  the  boundary  loads  at  ev¬ 
ery  observation  point  i  in  an  IPL  test  series  for  a  given  material  can  be  computed  by  numerical  integration 
since  all  tractions/  and  applied  boundary  displacements  u  required  for  the  computations  were  previously 
measured  and  stored.  Thus,  R 1  is  given  by 


where  i  =  On  the  other  hand,  the  energy  absorbed  by  the  specimen  is  by  linear  additivity 

=  Jvo/<p  (§(*))  4*  =  Cjjvol  Xj  (§'  (c,x))dx. 


(54) 


(55) 
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where  the  volume  integration  is  to  be  taken  in  a  Stieltjes  sense  over  the  plies.  When  doing  ply-level  analyses 
of  angle-ply  materials,  we  have  used  the  representation  given  by  Eq.  (28),  and  the  Einstein  summation  con¬ 
vention  is  in  effect.  But  <J>*  must  equal  by  energy  conservation,  therefore 

=  Cj\  X,(§i(c,x))dx.  (56) 

JJvol  1 

The  quantity  g'  (c,  x)  that  appears  in  Eq.  (55)  is  the  specimen  strain  field  corresponding  to  the  Ith 
observation  and  is  not  known  explicitly;  however,  it  must  satisfy  the  constitutive  relation  as  given  by  Eq. 
(49)  as  well  as  the  equilibrium  and  strain-displacement  equations,  i.e. 


diva'  =  0, 
and 


(57) 


§'  =  sym(gradu')-  (58) 

The  appropriate  boundary  conditions  are  free  surface  conditions  (no  tractions)  except  on  the  gripped 
surfaces  where  displacements  are  prescribed.  Equations  (54)  through  (58)  form  a  coupled  system  of  nonlin¬ 
ear  equations  whose  solution  places  restrictions  on  the  quantities  (c,  g,  O,  u)  . 

5.6  Method  of  Solution 

The  nature  of  the  system  of  Eqs.  (54)  through  (58)  suggests  solution  by  the  following  iterative  pro¬ 
cedure: 

1 .  An  initial  estimate  of  the  strain  field,  i.e.,  g'°  ( Jc)  is  obtained  by  solving  Eqs.  (57)  and  (58)  using 
the  constitutive  relation 


9  =  [C°]g  (59) 

where  [  C^]  is  the  stiffness  matrix  for  the  unstressed  material  and  can  either  be  obtained  from  the  literature 
(when  available)  or  estimated  from  mixture  theory.  The  problem  defined  here  is  linear  and  is  readily  solved 
by  finite  element  methods,  or,  in  certain  cases,  by  closed-form  procedures  (Muskhelishvili  (1953),  Mast 
(1973a,73b)). 

2.  Using  g'°  (x)  as  the  estimate  for  g'  (a,  x)  ,  an  estimate  of  the  parameter  vector  c  is  obtained  from 
Eq.  (56)  subject  to  the  constraints  that  each  component  of  the  imparted  energy  be  monotonically  increasing 
along  every  loading  path.  This  involves  solving  a  linear  or  quadratic  programming  problem  as  explained  in 
Section  5.7. 

3.  Equations  (57)^  (58),  and  (49)  are  solved  numerically  with  the  value  of  c  in  the  representation  for 
<|>  (Eq.  (28))  taken  as  c  ;  the  new  estimate  of  g  is  utilized  in  Eq.  (56)  to  obtain  an  updated  estimate  of  c; 
and  the  process  is  repeated  until  successive  estimates  differ  by  (hopefully)  a  sufficiently  small  amount. 

Note  that  except  for  the  initial  step,  determining  the  estimate  of  g'  (c,  x)  always  involves  solving  a 
set  of  nonlinear  boundary  value  problems.  Although  this  is  in  general  difficult  and  time  consuming,  as  a 
practical  matter,  unless  one  is  interested  in  venturing  deeply  into  the  nonlinear  regime  where  significant 
strain  redistribution  occurs  due  to  spatially  widespread  material  stiffness  changes,  sufficient  accuracy  should 
be  attainable  by  implementing  only  a  single  cycle  of  the  iteration  procedure.  Thus,  we  usually  need  only  to 
solve  the  linear  problem  described  in  Item  1  above,  and  the  difficulties  mentioned  above  should  be  of  little 
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concern.  If  the  need  does  arise  to  implement  more  than  one  cycle  of  the  iteration  procedure,  techniques  such 
as  homotopic  continuation  methods  (Watson  (1986))  applied  along  each  loading  path  are  available  to  help 
alleviate  the  difficulties. 

5.7  Initial  Estimates  of  the  Parameter  Vectors 

As  mentioned  previously,  g'°  (x)  is  computed  by  solving  a  linear  boundary  value  problem.  Indeed, 
using  the  fact  that  any  boundary  displacement  u  produced  by  the  IPL  can  be  represented  by  a  linear  com¬ 
bination  of  opening/closing  (uj),  sliding  (uq),  and  rotating  (uo)  displacements,  only  three  finite  element  anal¬ 
yses  need  be  performed  in  this  step  since  any  strain  field  £l0(x)  can  be  represented  as 

e'°  =  ^Te'0*  (60) 

k 

by  linear  superposition,  where  e'°*  (x)  are  the  strain  fields  corresponding  to  unit  boundary  displacements 

th  " 

in  the  k  direction. 

With  e'°  (x)  known,  the  quantity  JvoiXj  (g'0  (*))  in  Eq.  (56)  may  be  computed. 

Thus,  setting 

(61) 

Eq.  (56)  may  be  written  as  the  linear  set  of  equations 

[F]z=W,  (62) 

where  W=  (R\  ...,  Rn )  T;  zsc°,  and 


as  long  as  rank  ( [F] )  <  rank  (  [F] ) ,  where  [Fj  is  the  augmented  matrix 


(64) 


Equation  (62)  is  an  overdetermined  system  and  there  exists  no  value  of  z  that  satisfies  it  exactly.  Rath¬ 
er,  one  must  seek  to  minimize  the  norm  of  the  error  vector  e  =  (  W  -  [F]  z)  to  obtain  a  best  approxima¬ 
tion  to  the  solution.  This  is  the  usual  state  of  affairs  in  parameter  estimation;  indeed,  the  more  highly 
overdetermined  the  system,  the  better  the  situation  since  one  generally  wants  to  work  with  as  many  inde¬ 
pendent  observations  as  possible  to  filter  out  any  undesirable  stochastic  variations^ noise).  Only  if  enough 
observations  are  linearly  dependent  can  it  happen  that  rank  ( [FJ )  =  rank  ( [FJ )  <  m,  in  which  case 
there  exists  n  -  rank  (  [  FJ  )  solutions  to  Eq.  (62).  This  situation  should  not  arise  if  the  experimental  testing 
program  is  well  thought  out. 

Assuming  for  the  moment  that  Eq.  (62)(62)  is  indeed  overdetermined,  we  seek  to  minimize  the  norm 
of  e  =  W  —  [F]  z  subject  to  the  constraints  that 
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•  the  values  of  9  are  positive  at  the  node  points  of  the  interpolation  mesh,  i.e., 

z'  >  0,  (65) 

to  maintain  positivity, 

•  the  dissipated  energy  is  monotonically  increasing  along  each  radial  loading  path  in  displacement 
space,  i.e.. 


row  ( [FI),  -z>0 

*r  ~ 


(row  ( [F] ) f  +  j  -  row  ( [F] ) f)  ■  z  Z  0, 1  =  ip,  1 

where  Eq.  (66)  applies  to  each  loading  path,  ip  is  a  point  at  the  beginning  of  the  path,  and  iq  is  the  end 
point. 

Note  that  the  above  constraints  are  linear  inequality  constraints  in  z,  and  as  such  are  readily  dealt 
with.  Two  widely  used  measures  of  vector  magnitude  are  the  L2  (Euclidean)  and  L°°(Chebychev)  norms  giv¬ 
en  by 


WxWl*  =  (x]  +  x\+ ...+xl)2,  (67) 

and 

11*11  L-  =  max*  (jjf,|,  \x2\, ...,  |jcb|)  (68) 

respectively.  Minimizing  ||  «||  in  Eq.  (62)  under  the  l?  norm  subject  to  the  constraints  of  Eqs.  (65)  and  (66) 
is  a  problem  in  quadratic  programming  since  the  objective  function  ||  W  —  [F]  £||  2  is  quadratic  in  the  vari¬ 
ables  (Zp  Z2»  •••,  Z„) .  This  is  a  well-established  discipline  and  many  efficient  procedures  are  available  to 
deal  with  such  problems.  Minimizing  \e\  under  the  L°°  norm  is  just  as  straightforward,  and  perhaps  easier, 
because  there  exists  a  well-known  technique  due  to  Lagrange  for  converting  such  problems  into  linear  pro¬ 
gramming  problems  that  can  be  solved  using  the  simplex  method  or  one  of  its  variants.  A  highly  readable 
discussion  of  these  optimization  issues  is  presented  by  Gill  Murray,  and  Wright,  (1981). 
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Appendix  A 

Materials  Data  Base  Description 

ITEM 

RESIN  FIBER  LAYUP  SUPPLIER 

FILE 

Material  System  AS 1/3501 -6  (Thermoset) 

001 

3501-6  AS1  +/-  15  Hercules 

mt-1 

002 

+/-  30 

mt-2 

003 

+/-  60 

mt-3 

004 

+1-15 

mt-4 

Material  System  AS4/3501-6  (Thermoset) 

005 

AS4  +/-15 

mt-5 

006 

+/-30 

mt-6 

007 

+/-  60 

mt-7 

008 

+1-15 

mt-8 

Material  System  AS6/3501-6  (Thermoset) 

009 

AS6  +/-15 

mt-9 

010 

+/-30 

mt-10 

Oil 

+/-60 

mt-11 

012 

+1-15 

mt-1 2 

Material  System  IM 6/3 50 1-6  (Thermoset) 

013 

IM6  +/-I5 

mt-1 3 

014 

+/-30 

mt-14 

015 

+/-60 

mt-15 

016 

+1-15 

mt-16 

Material  System  IM7/3501-6  (Thermoset) 

017 

IM7  +/-15 

mt-1 7 

018 

+/-30 

mt-1 8 

019 

+/-60 

mt-19 

020 

+1-15 

mt-20 

PROJECT 

NAV-SEA 
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ITEM 

RESIN  FIBER  LAYUP 

SUPPLIER  FILE 

Material  System  HMS4/3501-6  (Thermoset) 

021 

3501-6  HMS4  +/-15 

Hercules  mt-21 

022 

+/-30 

mt-22 

023 

+1-60 

mt-23 

024 

+1-15 

mt-24 

Material  System  HS  12000/5245  (Thermoset) 

025 

5245  HS  12000  +/-15 

Celanese  mt-25 

026 

+1-15 

mt-26 

Material  System  HS  12000/5208  (Thermoset) 

027 

5208  +/-15 

mt-27 

028 

+1-15 

mt-28 

Material  System  6000/5208  (Thermoset) 

029 

CELION  6000  +/-  15 

mt-29 

030 

+/-15 

mt-30 

Material  System  AS4W/3501-6  (Thermoset) 

031 

3501-6  AS4W  0,90,0%z 

FMI  mt-31 

032 

0,90,5%z 

mt-32 

033 

0,90, 13%z 

mt-33 

Material  System  IM6/PEEK  (Thermoplastic) 

034 

PEEK  IM6  +/-15 

ICI  mt-34 

035 

+/-30 

mt-35 

036 

+/-60 

mt-36 

037 

+/-15 

mt-37 

Material  System  AS4/PEEK  (Thermoplastic) 

038 

AS4  +/-15 

mt-38 

039 

+/-30 

mt-39 

040 

+/-  60 

mt-40 

041 

+/-15 

mt-41 

PROJECT 
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ITEM  RESIN  FIBER  LAYUP  SUPPLIER  FILE 
Material  System  E-glass/Epoxy  (Thermoset) 


042 

Epoxy  E-glass  Cloth  Fiberite 

mt-42 

043 

Magnaweave 

mt-43 

044 

. 

mt-44 

Material  System  B238/Polyester  (Thermoset) 

045 

Polyester  B238  0,90  Proform  Fab. 

mt-45 

046 

.  ((0,90)(+/-45))2 

mt-46 

047 

.  (0,90)(+/-45)2(0,90)  . 

mt-47 

048 

(+/-45)4 

mt-48 

Material  System  CDB340/Polyester  (Thermoset) 

049 

CDB340  (0,+/-45) 

mt-49 

050 

(+/-45.90) 

mt-50 

Material  System  DB170/Polyester  (Thermoset) 

051 

DB170  (0,90)5 

mt-51 

052 

(+/-45,0,90)2 

mt-52 

053 

.  (0,90)(+/-45)3(0,90)  . 

mt-53 

054 

(+/-45)s 

mt-54 

Material  System  Quartz/Polybutadiene  (Thermoplastic) 

055  Polybutadiene  Quartz  +/-0  Brunswick 

mt-55 

056 

+/-15 

mt-56 

057 

+/-30 

mt-57 

058 

+/-  45 

mt-58 

059 

+/-60 

mt-59 

060 

+1-15 

mt-60 

061 

+/-90 

mt-61 

PROJECT 


RAS 
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ITEM 

RESIN  FIBER  LAYUP  SUPPLIER 

FILE 

PROJECT 

Material  System  Quartz/Poly  butadiene  5%  Powder  Carbon  (Thermoplastic) 

062Polybutadiene,5%C,  Quartz  +/-0  Brunswick 

mt-62 

• 

063 

+/-15 

mt-63 

064 

+/-30 

mt-64 

065 

+/-45 

mt-65 

066 

+/-60 

mt-66 

067 

+/- 75 

mt-100 

068 

+/-90  .  mt-101 

Material  System  IM7/HBRF-55A  (Thermoset) 

069 

HBRF-55A  IM7  +1-22  Hercules 

mt-70 

SP 

070 

+/-35 

mt-71 

071 

+/-55 

mt-72 

072 

+/-68 

mt-73 

073 

(+/- 22,90) 

mt-74 

074 

(+/-68.0)  .  mt-75 

Material  System  IM7/HBRF-55A  (Thermoset) 

075 

HBRF-55A  IM7  +1-22 

mt-102 

076 

.  (-55,+35,-35,+35) 

mt-103 

077 

.  (+35, -55, +55, -55) 

mt-104 

078 

+/-68 

mt-105 

079 

(90.+/-22) 

mt-106 

080 

(0.+/-68) 

mt-107 

Material  System  Kevlar/Polyether-Sulfone  (Thermoplastic) 
081Polyether-SulfoneKevlar  49  Fabric  ICI  mt-108  Chem. 

Material  System  TORAY  300/Polyether-SuIfone  (Thermoplastic) 

082  .  TORAY  300  .  .  mt-109 

Material  System  E-glass/Polyether-Sulfone  (Thermoplastic) 

083  .  E-glass  .  .  mt-110 

Material  System  Graphite  A004/Epikot  828  (Thermoset) 

084  Epikote  828  A004  (0,90,90/0)  TTCP  mt-111  TTCP 

Material  System  S-2  glass/Shell  828(Thermoset) 

085  Shell  828  S-2  glass  (0,90)5  mt-112 
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ITEM 

RESIN 

FIBER 

LAYUP 

SUPPLIER 

FILE 

PROJECT 

Material  System  RYTON  PPS/GRAPHITE  (AC  40-66)  (Thermoset) 

086 

Ryton  PPS 

graphite 

+/-15 

Phillips 

mt-113 

Boeing 

087 

• 

• 

+/-30 

. 

mt-114 

• 

088 

• 

+/-60 

. 

mt-115 

. 

089 

♦ 

. 

+/-75 

. 

mt-1 16 

. 

Material  System  IM7/HBRF-55A  (Thermoset) 

090 

HBRF-55A 

IM7 

+/-  22 

Hercules 

mt-1 17 

SP-D5 

091 

. 

. 

+/-35 

mt-1 18 

. 

092 

. 

. 

+/-  55 

mt-1 19 

093 

• 

+/-68 

mt-1 20 

094 

• 

• 

(+/- 22,90) 

mt-121 

• 

095 

. 

. 

(+/-68.0) 

mt-1 22 

Material  System  AS4/PEEK  (APC-2)  (Thermoplastic) 

096 

PEEK 

AS4 

+/-15 

• 

mt-1 23 

097 

. 

+/-30 

. 

mt-1 24 

098 

. 

+/-45 

. 

mt-1 25 

099 

. 

+1-60 

. 

mt-1 26 

100 

. 

+1-15 

• 

mt-1 27 

101 

(0,90) 

. 

mt-1 28 
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Appendix  B 

The  In-Plane  Loader  System 


The  IPLS,  shown  schematically  in  Fig.  61,  is  an  automated  testing  facility  whose  primary  function  is 
to  produce  reliable  data  from  which  the  energy  dissipation  function  can  be  deduced  with  a  minimum  of  an¬ 
cillary  hypotheses. 


R.  LOADING  I 


B.  HVDRRULIC  RCTURTORS 

C.  MOURBLE  HERD 

D.  FIXED  HERD 

E.  SPECIMEN  GRIPS 
F-  SPECIMEN  LOADER/UNLORDER 

G.  TEST  SPECIMEN 

H.  OIDEO  CRMERR 

I.  UIDEO  MONITOR  FOR  POSITIONING  DEVICE 

J.  MRIN  TERMINRL 
X.  COMPUTER 

L.  GRRPHICS  PROCESSOR 

M.  GRRPHICS  TERMINRL 


FRRME 


Fig.  B 1  —  Schematic  representation  of  the  In-Plane  Loader  System 


B.l  System  Components 
B.1.1  The  In-Plane  Loader 

This  is  a  sophisticated  automated  testing  machine  capable  of  subjecting  test  specimens  to  arbitrary 
displacement-controlled  in-plane  loading  corresponding  to  combinations  of  pulling/pushing,  sliding,  and 
rotating  boundary  displacements. 

B.l. 2  The  Computer  System 

The  computer  serves  as  a  controller  by  sending  commands  to  the  hydraulic  actuators  of  the  IPL 
through  digital-to-analog  converters,  and  by  receiving  gauge  measurements  from  the  IPL  via  the  analog-to- 
digital  converters.  It  is  also  used  for  real-time  data  collection  and  reduction  during  the  tests  as  well  as  for 
various  post-testing  analyses.  Among  the  latter  are  the  determination  of  the  energy  dissipation  function  and 
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consequent  constitutive  behavior  for  the  material  and  finite  element  analysis  of  specific  structural  compo¬ 
nents. 

B.1.3  The  Graphics  Processor 

The  graphics  processor  is  used  to  facilitate  the  pre-  and  post-processing  requirements  of  the  finite  el¬ 
ement  analyses.  Specifically,  it  is  used  in  conjunction  with  solid  modeling  and  rendering  software  to  produce 
finite  element  models  of  the  specimen  itself  and  specific  structural  components.  It  is  also  used  to  generate 
rendered  images  showing  the  spatial  distribution  of  dissipated  energy. 

B.2  In-plane  Loader  Description 

Although  the  IPL  has  been  described  previously  (Mast  et  al.  (1983)),  its  design  has  since  evolved  suf¬ 
ficiently  to  warrant  an  updated  description.  The  in-plane  loader  in  its  current  incarnation  is  shown  in  Fig. 

B 1  with  a  more  detailed  portrayal  in  the  vicinity  of  the  specimen  grips  shown  in  Fig.  B2.  Basically,  it  con- 


Fig.  B2  —  In-Plane  Loader’s  closeup  view  of  the  grip  area 


sists  of  three  independent  computer  controlled  hydraulic  actuators  connected  to  a  movable  head.  The  load¬ 
ing  path  of  each  actuator  can  be  independently  prescribed,  thereby  producing  arbitrary  combinations  of 
opening/closing,  sliding,  and  in-plane  rotating  displacements. 

Each  actuator  unit  is  composed  of 

•  2-Kip  hydraulic  linear  actuator  with  a  6-in.  stroke, 

•  DC  powered  linear  displacement  transducer  (DCLVDT)  which  is  mounted  on  the  actuator  and  mea¬ 
sures  ram  motion  relative  to  the  actuator  cylinder, 

•  operational  amplifier  that  computes  the  difference  (the  displacement  error)  between  the  actual  ram 
displacement  and  the  programmed  one, 
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•  electrically  controlled  hydraulic  servo  valve  that  receives  the  displacement  error  output  of  the  opera¬ 
tional  amplifier  and  delivers  controlled  hydraulic  power  to  the  actuator,  forcing  the  ram  to  displace 
the  amount  necessary  to  minimize  the  displacement  error, 

•  2- Kip  strain  gauge  load  cell  mounted  on  the  ram  end,  and 

•  instrument  amplifier  to  amplify  the  millivolt  output  of  the  load  cell. 

Relative  displacements  between  the  two  grips  are  monitored  by  a  specially  designed  six-component 
DCLVDT  mounted  on  the  specimen  side  of  the  movable  head.  A  three-component  strain-gauged  load  cell 
is  used  to  monitor  the  movable  head  reaction  forces. 


Fig.  B3  —  Schematic  of  the  equivalence  between  the  (a)  actuator  motions  and  (b)  the  three  modes  of  motion .  as 
well  as  (c)  the  process  of  composing  an  arbitrary  loading  combination,  from  linear  combinations  of  the  three 
basic  motions:  (d)  shear,  (e)  opening/closing  and  (f)  rotation 


Two  actuator  configurations  deserve  special  attention.  The  first  is  when  any  of  the  actuators  is  parallel 
with  another,  while  the  second  occurs  when  the  lines  of  action  of  the  actuators  intersect  at  a  common  point. 
Both  of  these  configurations  prevent  the  application  of  the  full  range  of  opening/closing,  sliding,  and  rota¬ 
tion,  and,  consequently,  were  avoided  in  the  design. 

The  hydraulic  grips  are  computer  controlled  and  operate  like  a  pair  of  pliers  whose  handles  are  forced 
together  by  a  compact  servo-controlled  jack  built  into  the  head.  Specimens  are  inserted  in  the  grips  from  a 
loading  magazine  via  a  mechanical  feeding  device.  This  ensures  uniformity  in  the  gripping  conditions  and 
allows  for  a  rapid  throughput  (currently  340  specimens/hour).  A  hydraulic  power  control  unit  consisting  of 
an  operational  amplifier,  a  pressure  transducer,  and  an  electrohydraulic  valve,  maintains  a  constant  gripping 
pressure  as  the  test  proceeds.  The  magnitude  of  the  gripping  force  is  selected  to  be  sufficient  to  prevent  spec¬ 
imen  slippage,  and  yet  not  so  great  as  to  cause  failure  at  the  grips. 

The  digital  imaging  system  is  provided  to  give  the  operator  an  efficient  and  accurate  method  of  per¬ 
forming  initial  specimen  geometry  measurements  and  calibrations.  A  video  camera  with  a  long  focal  length, 
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positioned  by  a  motor  driven  three-axis  stage,  is  mounted  over  the  specimen  and  provides  a  magnified  (15x) 
view  of  the  specimen  on  a  TV  monitor.  Horizontal  and  vertical  cursor  lines  are  superimposed  on  the  TV 
image  and  can  be  used  to  mark  various  features  such  as  the  initial  position  of  the  notch  tip.  The  cursor  digital 
output  signals  are  continuously  monitored  and  saved  on  command. 

The  block  diagram  of  the  entire  system  is  shown  in  Fig.  B4  together  with  a  schematic  representation 
of  the  information  paths  from  and  to  the  IPL  relative  to  the  computer  controller.  Twenty-one  measuring 
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Fig.  B4  —  Block  diagram  of  the  In-Plane  Loader  System 


channels  are  interfaced  with  a  multiplexer  connected  to  the  analog-to-digital  converter  which  is  connected 
to  the  computer  through  DMA  (direct  memory  access)  via  the  computer  bus.  Similarly,  the  digital-to-analog 
converter  has  its  output  connected  to  the  servo  valves  of  the  actuators  and  its  input  communicates  with  the 
computer  via  the  bus  in  a  DMA  fashion. 
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Appendix  C 

A  Vector  Space  Representation  of  the  Optimization  Scheme 


A  vector  space  geometrical  representation  of  the  optimization  scheme  provides  for  a  simple  explana¬ 
tion  of  the  optimization  procedure.  Assume  for  the  sake  of  simplicity  (and  because  we  can’t  represent  enti¬ 
ties  having  dimensions  higher  than  three  on  paper)  that  the  quantity  [X]  c  can  be  represented  as  a  vector 
c  in  a  two-dimensional  space  X  with  Xj  and  X2  being  the  corresponding  basis  vectors  as  shown  in  Fig. 

C 1 .  The  columns  of  the  array  [X]  contain  the  components  of  the  basis  vectors  Xj  and  X2  in  terms  of  a  higher 
dimensional  space  D,  the  higher  dimensions  of  which  are  denoted  by  the  axis  N  normal  to  the  space  X. 

Assume  now  that  another  vector  d  lies  in  the  higher  dimensional  space  D  of  which  the  space  X  is  a 
subspace  as  shown  in  Fig.  15.  The  difference  of  these  two  vectors  then  lies  in  this  higher  dimensional  space 
and  is  another  vector  denoted  as  e  which  expresses  the  error  between  the  two  vectors,  so  that 

e  =  d-c  .  (Cl) 

This  can  be  rewritten  in  the  matrix  form 

[X]c  +  g  =  d.  (C2) 

The  objective  of  the  process  is  to  determine  a  vector  CQ  such  as  the  magnitude  (or  norm)  of  the  vector 
e  takes  its  minimum  value  according  to 


min  ||  e||  =  eQ.  (C3) 

In  the  case  of  the  usual  Euclidean  norm,  i.e.,  the  square  root  of  the  sum  of  the  squares  of  the  components, 
Cq  turns  out  to  be  the  projection  of  the  vector  d  of  space  D  onto  the  space  X,  and  therefore  e0  is  normal  to 
this  space. 

Fig.  Cl(b)  deals  with  the  effects  of  inequality  constraints  case.  Suppose  Mx ,  M2  are  the  vectors  cor¬ 
responding  to  the  matrices  M  j  and  M 2  ,  which  also  lie  in  the  space  X.  Then,  the  inequality  constraint  of 
the  form 


[Mleso  (C4) 

specifies  that  the  only  acceptable  solutions  for  c  must  lie  in  the  region  of  X  (shaded  region)  defined  by  the 
interior  of  the  two  lines  which  also  lie  in  X  and  are  normal  to  the  vectors  A/j ,  M2  respectively. 

However,  as  can  be  seen  in  Fig.  Cl(c),  the  vector  c0,  which  exactly  satisfies  Eq.  (22)  and  minimizes 
e  to  a  value  eQ,  may  lie  outside  the  region  specified  by  the  constraints.  In  this  case  the  quadratic  program¬ 
ming  methodology  determines  a  vector  c'Q,  which  satisfies  the  constraints  by  being  the  projection  of  cQ 
onto  the  closest  boundary  of  the  shaded  region.  The  magnitudes  of  the  components  of  c'Q  in  the  space  X 
will  then  be  the  desired  coefficients  c,.  Finally,  this  process  is  valid  when  the  spaces  X  and  D  have  arbitrary 
dimensions  so  long  as  Dim  (X)  <  Dim  (D)  and  so  obviously  applies  in  our  case. 
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Fig.  C 1  —  Geometrical  representation  of  the  optimization  procedure  describing  ( a )  the 
solution  for  Eq.  (C2),  ( b )  for  the  constraints ,  and  ( c )  for  both  equation  and  constraints 
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